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ABSTRACT 


A  preprocessor  is  designed  to  extract  a  set  of  features 
that  enhance  natural  clustering  by  removing  extraneous 
information.  The  design  removes  time  shift  and  scale 
dependence  by  taking  advantage  of  invariant  properties  of  a 
Fourier  transform  followed  by  a  mellin  transform.  The 
preprocessor  is  realized  using  an  FFT  and  a  Mellln 
transform  with  a  conventional  error  correction  term.  The 
error  term  proves  to  be  indeterminate,  but  the  error's 
bound  is  identified  as  the  envelope  for  Pellin  correction 
terms.  Properties  of  the  Mellin  transform  are  employed  to 
modify  the  signal  so  that  the  error  correcting  is  no  longer 
required.  The  resulting  algorithms  are  tested  with 
variously  scaled  inputs  for  which  closed  form  solutions  are 
known.  With  a  verified  modification  in  place,  the 
preprocessor  produces  features  that  are  invariant  to 
shifting  and  scaling,  while  retaining  enough  information  to 
classify  canonic  shapes.  A  method  of  improving  performance 
is  introduced. 
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I.  INTRODUCTION  AND  BACKGROUND 


Pattern  classification  is  the  assignment  of  a  physical 
object  or  event  to  one  of  several  prespecified  catagories 
and  is  the  result  of  an  incomplete  theory  of  perception. 
Although  many  transducers  are  available  for  converting 
light,  sound,  temperature,  reflected  radar  signals,  etc., 
to  electrical  signals,  the  ability  of  machines  to  perceive 
cr  to  recognize  their  environment  remains  very  limited.  In 
the  structured  world  of  communcation s  engineering,  signals 
are  designed  to  be  detectable  and  differentiable.  A  much 
more  difficult  problem  presents  itself  when  sensing  an 
environment  through  a  transducer  and  recognizing  or  even 
classifying  the  elements  of  that  environment  on  the  sensed 
characteristics  of  the  transducer's  electrical  output. 
Pattern  recognition  can  be  considered  a  complex 
communications  problem  (for  example,  attempting  to  teach  a 
machine  to  decode  signals  encoded  by  nature).  It  is 
possible  to  alter  the  transducer's  output  to  facilitate 
object  classification,  but  determining  how  to  alter  that 
output  is  not  a  simple  task.  The  main  elememts  of  a 
classification  system  is  shown  in  figure  1  [1].  The 
transducer  senses,  actively  or  passively,  a  set  of 
characteristics  belonging  to  the  object.  These 
characteristics  can  never  be  a  complete  description  of  an 
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A  Classification  System 
Figure  1 


object,  but  represent,  hopefully,  enough  Information  to 
classify  the  object  as  belonging  to  one  of  a  number  of 
classes.  For  Instance,  temperature  Is  a  characteristic  of 
the  object  and  a  class,  but  this  feature  Is  of  little  value 
unless  It  differs  in  some  way  from  objects  belonging  to 
other  classes,  while  the  set  must  Include  characteristics 
that  are  common  among  that  class.  The  preprocessor  (or 
feature  extractor)  aims  to  reduce  the  data  by  measuring  or 
quantifying  certain  properties  that  distinguish  the  sensed 
object  as  belonging  to  one  class  and  not  to  others.  This 
can  be  done  by  discerning  key  features  or  using  the  imput 
to  generate  another  set  of  features  optimized  by  some  rule. 
The  values  of  each  of  these  features  is  then  passed  to  the 
classifier,  which  evaluates  these  features  to  assign  the 
object  to  a  class. 

With  varied  success,  machine  pattern  classification  has 
been  applied  to  a  large  range  of  problems/disciplines. 
Fields  where  it  is  particularly  common  Include  optical 
imagery,  acoustic  signal  processing,  radiology,  radio 
astronomy,  and  electronic  warfare,  to  name  a  few.  Work  in 
many  of  these  fields  was  reviewed  during  the  development  of 
this  thesis  and  the  results  derived  and  demonstrated  here 
may,  in  turn,  be  applicable  to  the  field  in  general.  This 
effort  has  been  directed  toward  designing  a  preprocessor  to 
produce  a  set  of  features  that  are  invariant  to  information 
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known  to  be  superfluous  to  classification,  but  that  retain 
enough  Information  to  classify  an  object.  The  object  has 
been  sensed  by  a  transducer  and  has  been  represented  as  an 
empirically  derived,  univariate  time  series.  Such  a  series 
would  be  the  form  of  data  available  from  range  only  radar 
return  which  is  specifically  what  the  preprocessor  is 
designed  to  handle.  Returning  to  figure  1,  the  object  has 
an  infinite  set  of  characteristics  (here  portrayed  as  an 
infinite  series  of  discrete  values).  The 
transducer/receiver  has  collected  some  characterist ics  of 
the  target  objective  in  the  presence  of  noise.  This 
information  is  relayed  as  a  set  of  discrete  samples  (hi) 
from  a  band  limited  signal.  The  preprocessor  is  designed  to 
determine  and  code  revelent  information  (Hj)  for  the 
classifier.  If  this  task  was  done  well,  classification 
becomes  a  trivial  problem.  On  the  other  hand,  if  the 
classifier  becomes  ideal  (capable  of  resolving  an  infinite 
number  of  characteristics  in  noise)  the  preprocessor  design 
begins  to  look  like  a  wire.  The  distinction  between  the 
preprocessor  and  the  classifier  is  arbitrary  from  an 
analytical  point  of  view.  When  designing  a  classification 
system  functionally,  a  difference  is  enforced.  The 
classifier  has  little  concern  for  how  the  features  are 
developed,  but  seeks  to  efficiently  use  those  provided  to 
guess  the  class  of  the  target  object.  The  preprocessor  is 
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problem  dependent,  needing  to  produce  an  optiiral  set  of 
features,  Ej ,  from  the  sensed  data  hi. 

A.  FEATURE  EXTRACTION 

For  the  purposes  of  this  paper,  two  generic  approaches 
to  feature  extraction  are  defined.  The  first,  a 
classification  approach,  was  described  above.  The  second,  a 
descriptive  approach,  tries  to  define  the  object  in  terms 
of  the  objects'  structural  features.  This  system  might 
recognize  a  car,  for  example,  by  breaking  up  the  visual 
picture  into  canonic  shapes,  and  comparing  this  to 
previously  specified  canonic  class  models.  The  perceived 
structure  of  the  physical  object  is  maintained  and  should 
reflect  the  structure  of  the  object  itself.  This  approach 
could  be  robust  to  temporary  changes  in  the  object  itself. 
In  the  car  example,  knowing  that  at  one  end  of  the  car  the 
trunk  can  be  opened  or  closed  allows  the  device  to  take 
this  factor  into  account.  Another  important  advantage  to 
descriptive  techniques  is  that  the  class  characteristics 
may  be  entered  or  specified  without  collecting  actual 
transducer  generated  data  to  train  the  machine. 
Unfortunately,  the  problem  of  designing  a  machine  to 
analyse  a  visual  scene  to  produce  a  structural  description 
has  proved  to  be  quite  difficult.  Object  description  from  a 
univariate  time  series  is  even  more  difficult,  and  if  the 
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radiation  sensed  by  the  transducer  is  not  froir  the  visual 
spectrum,  the  task  rapidly  approaches  the  impossible.  For 
these  reasons  the  approach  taken  was  the  classification 
approach  (to  reduce  the  signal  to  a  set  of  orthogonal 
features  that  do  not  uniquely  reflect  the  structure  of  the 
object,  but  do  retain  sufficient  information  to  classify 
the  object ) . 

This  paper  excludes  a  detailed  description  of  the 
transducer  specification.  The  problem  of  the  classifier 
itself  is  viewed  as  one  of  partitioning  the  feature  space 
(Hj)  into  regions;  one  region  for  each  category.  Ideally, 
this  partitioning  should  be  arranged  so  that  none  of  the 
decisions  are  ever  wrong.  When  this  cannot  be  realized,  at 
least  the  probability  of  error  should  be  minimized  or  the 
average  cost  of  errors  minimized.  The  problem  is  one  within 
statistical  decision  theory,  inowledge  of  the  object 
classes  (the  transducer  and  the  classifier)  are  required  to 
design  the  preprocessor,  which  is  the  topic  of  this  thesis. 
The  preprocessor  designed  and  built  here  generates  features 
from  a  range  only  radar  video  signal.  These  features  are 
used  by  a  general  Bayesian  learning  classifier.  The 
supervised  learning  general  Bayesian  classifier  approaches 
the  problem  by  taking  a  series  of  incoming  sets  of  features 
labeled  as  to  their  class.  From  the  data,  an  a  posteriori 
density  is  computed.  Each  successive  set  of  training  data 
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is  used  to  refine  the  densities'  statistics,  rfhen  the 
classifier  has  been  trained  on  N  classes,  the  features  are 
irodified  to  separate  the  class  volurres  in  an  optimal  way 
and  to  reduce  the  number  of  features  to  one  less  than  the 
number  of  possible  classes.  The  feature  vectors  of  class 
members  are  clustered  about  a  simplex  point  and  likely 
boundaries  are  set  up  allowing  classification  of  the  object 
as  belonging  to  one  of  the  classes,  or  of  an  unkown  class. 
An  A'  simplex  is  a  collection  of  N  points  in  (N-l)  space 
where  the  distance  between  any  two  of  the  points  is  equal 
to  the  distance  between  any  other  two.  Thus  a  three  class 
problem  transformed  into  a  three  simplex  in  a  two 
dimensional  plane  produces  clustering  of  the  three  class's 
about  the  vertices  of  an  equilateral  triangle.  The  simplex 
coordinates  are  the  reduced  feature  vectors,  generalized 
from  the  training  data.  In  a  controlled,  simple  problem  the 
classifier  works  well,  but  when  encountering  real  problems 
one  class's  feature  space  will  intertwine  another's,  making 
it  much  more  difficult  to  obtain  separation  in  a  meaningful 
way.  The  goal  of  the  preprocessor  is  to  present  key 
features  that  determine  class  for  subsequent  optimization 
by  the  classifier. 
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B.  JOUR  HR  -  MELLIN  PREPROCESSING 


Common  to  all  univariate,  time  series  classification 
problems  are  several  variables  that  interfere  with  the 
recognition  process.  Assuming  discrete  data  processing  is 
used,  these  are  addressed  in  the  following  order* 
windowing,  framing,  scaling,  sampling  rate,  quantization 
noise,  and  sufficient  Information.  For  real  processing,  the 
imput  waveform  is  not  sampled  for  all  time.  It  is  sampled 
for  a  period  of  time.  This  windowing  of  the  data  corrupts 
the  resulting  spectrum  in  two  ways  [2].  First  it  introduces 
a  periodicity  (the  inverse  of  the  window  length)  and 
resulting  aliasing  to  the  otherwise  infinite  spectrum,  and 
further  distorts  the  spectrum  by  a  convolution  with  the 
spectrum  of  the  window  itself.  3oth  of  these  effects  will 
color  all  of  the  data  in  the  same  way  and  so  can  be 
accounted  for  by  deterministic  methods.'  Framing  can  be 
considered  characteristic  of  poor  synchronization,  whereby 
the  pattern  of  concern  is  not  position  stable  with  repect 
to  the  window  as  shown  in  figure  2a.  It  seems  that  even  in 
human  optical  recognition,  the  eye  tends  to  center  the 
pattern  prior  to  recognition,  unless  trained  otherwise. 
Scaling  is  that  property  whereby  the  object  field  may  vary 
in  scale  or  aspect  angle,  in  one  dimension  or  several 
dimensions  as  in  figure  2b.  Before  the  analog  signal  is 
sampled,  it  must  be  fed  through  a  sharp  low  pass  filter. 
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Shifting  and  Scaling  Variation 
Tigure  2 


because  no  higher  frequency  noise  can  be  present  without 
being  folded  onto  the  valid  data.  Quantization  noise,  due 
to  the  requirement  to  round  off  each  sample  to  some 
discrete  level,  Is  treated  the  same  as  round  off  error  In 
numerical  processing  [3].  In  all  of  these  problems 
discussed  to  this  point,  the  effect  of  this  processing  is 
to  mask  the  actual  feature  vectors,  making  the 
classification  system  less  sensitive  to  valid  pattern 
characteristics.  In  all  recognition  problems,  it  is  assumed 
that  there  is  sufficient  information  present  for  a  pattern 
to  be  detected  and  identified  by  the  system.  This  means 
that  there  is  sufficient  variability  between  the  classes, 
but  sufficient  similarity  between  those  patterns  of  a 
single  class  to  classify  each  pattern  in  terms  of  those 
classes . 

A  verifiable  goal  of  this  preprocessor  is  to  produce  a 
set  of  features  that  are  invariant  to  shifting  and  scaling 
changes.  An  approach,  figure  3,  has  been  proposed  and  used 
[4-10]  .  The  preprocessor  consists  of  putting  the  sampled 
data  successively  through  two  transforms,  a  fast  Fourier 
transform  and  a  discrete  Mellin  transform.  Integral 
transforms,  the  Fourier  and  Mellin  included,  develop 
naturally  from  the  solution  of  simple  problems  in  potential 
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theory  [11-13].  The  Fourier  Integral  transform  of  the 
waveform  h(t), 

H(f)  =»  ^ ‘AM  K(f,  (i) 

where  I=exp  [-2ir*f t] ,  is  a  principle  analytical  tool  in  such 
diverse  fields  as  linear  systems,  optics,  probability 
theory,  quantum  physics,  and  signal  analysis  [13].  Its 
purpose  in  the  preprocessor  is  twofold,  hut  relies  on  a 
single  characteristic.  The  magnitude  of  the  Fourier 
transform  is  invariant  to  shifting,  h(t-a). 

I  Hal  I  =  /  H<t)l  (2) 

This  characteristic  removes  the  effects  of  framing 
inaccuracies,  and  also  permits  the  averaging  of  successive 
looks  or  pulses  of  data  to  improve  feature  resolution,  but 
removes  much  of  the  information  about  a  signal's  structure 
as  discussed  later.  A  discrete  Fourier  transform, 

4  )  *  O ,  t  j  T./  •»  »  ,  M  —  / 

HU)  -*  H  'VN  «  <S/  >!  2J  1  (3) 

a  €  (4) 

has  an  equivalent  identity,  but  it  is  only  exact  for  shifts 
of  integer  values. 
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(5) 
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f  3")  C  Ac**--*  )"\  J 

~  j  H(m)  *  |H  (** )  j 

Shifts  of  other  than  integer  values  result  in  errors  that 
depend  not  only  on  the  shift,  hut  on  qualities  of  the 
sampled  waveform  itself,  h(t). 

The  Mellin  transform  is  an  integral  transform  with  the 
kernel,  K. 

K  *  i  (6) 

H<->  *  j* l(t)  t'"  Jt  (?) 

is  the  Mellin  transform  with  respect  to  the  complex 
parameter  s=r-J24hm.  Several  simple  substitutions  relate 
this  to  more  common  analytical  tools.  Exponentially  warping 
t*exp[x]  changes  the  appearance  of  the  integral  to  what  is 
often  called  a  double  sided  Laplace  transform  [14]  , 

Hu)  =*  C  £ (e.*)  ^ 

The  transform  is  invariant  to  t  domain  scaling  when  taken 
with  respect  to  the  imaginary  part  alone, 

H a  f  “£  '<+)  a*  ( 9 ) 
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Equation  (S)  is  recognized  as  the  Fourier  intergral  of  an 
exponentially  distorted  waveform  h'(x).  The  modulus  of  the 
expression  on  the  right  is  the  magnitude  of  the  Fourier 
transform  of  the  exponentially  distorted  time  function.  The 
property  to  he  exploited  in  a  Fourier-Mellin  (FM) 
preprocessor  is  that  the  modulus  of  the  transform  in  s,  is 
invariant  to  t-scaling.  Given  a  time  waveform  h(t),  its 
Mellin  transform  H(s)  is  given  as  equation  (7).  Scaling  the 
entire  t-domain  hy  k, 

Letting  r=t/k,  and  remembering  that  s  is  imaginery, 
iJ  *  A*  H<s) 

IA'Vm  I  *  I  He sjI  (11) 

Much  effort  and  detail  is  spent  implementing  equation  (7) 
digitally  in  Chapters  II  and  III.  The  rest  of  this  chapter 
is  devoted  to  providing  some  required  background  on  the 
discrete  versions  of  the  Fourier  transform  and  some 
properties  upon  which  the  FM  preprocessor  depends.  The 
treatment  here  is  brief,  being  mainly  a  review  of  basic  FFT 
concepts  and  as  such  may  be  skipped  without  loss  of 
content. 

A  discrete  Fourier  transform  is  not  computationally 
efficient  and  so  leads  to  impractlcally  long  processing 
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times.  The  fast  Fourier  transform  (FFT)  efficiently 
computes  the  discrete  transform  and  is  used  in  the 
preprocessors  built  for  this  thesis.  Other  properities  of 
the  FFT  should  be  presented  before  getting  into  the 
detailed  design  of  the  preprocessor.  The  first  and  last  are 
concerned  with  symmetry.  If  h(m)  is  real,  as  in  the  case  of 
the  sampled  data,  then  the  frequency  spectrum  of  that  data 
is  even,  ! Ee(n) ! =i He (-n ) ! .  H(n)  has  a  real  part  and  an 
imaginary  part,  P.e(n)  and  Im(n),  while  h(m)  has  an  even 
he(m)®he (-m ) ,  and  an  odd  ho (m) =-ho (-m)  part. 

_  i  /  X'tt  fT\  /WS  | 

*«  C/"0  *  /  \  ^  ) 

m-i 

L~'C~)9  £  JL.  x;H  ( -i? /»  ^ ) 

'  N  )  (12) 

The  odd  part  of  h(n)  times  the  cosine  kernel  summation,  and 
the  even  part  of  h(n)  times  the  sine  kernel  summation  are 

both  zero.  H(n)  can  now  be  seen  to  have  an  even  and  odd 

part.  Talcing  the  magnitude  of  an  odd  function  makes  it 
even,  proving  that  the  Fourier  transform  of  a  real  series 
is  a  spectrum  whose  magnitude  is  symmetric  about  f=0.  This 
is  of  course  true  for  both  the  integral  and  discrete 
Fourier  transforms.  The  second  line  of  symmetry  is  an 
effect  of  discretizing  the  signal  and  its  spectrum  for  a 
discrete  transform  as  evinced  by  its  theoretical 

development.  Before  proceeding  though,  the  convolution 
theorem  is  required. 
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The  convolution  of  the  two  functions  h(t)  and  g(t)  is 
defined  as  the  familiar  integral, 
jft)  =  gU-i-'  Jr  - 

The  relationship  between  convolution  and  the  Fourier 
integral  is  very  important  to  modern  analysis  and 
contributes  to  making  the  Fourier  transform  a  key  analytic 
tool.  The  convolution  theorem  states  that  if  h(t)  has  a 
Fourier  transform  F(f)  and  g(t)  has  the  Fourier  transform 
G  ( f  ) ,  then  h(t)*g(t)  has  the  Fourier  transform  H(f)G(f). 

y  (t)  -  C«;  Ctj 


3-l^ctjJ  =  H</)  <3,(0 


(14) 


Proving  this,  the  Fourier  integral  is  used  directly  on  the 
convolution  integral. 

Y«)  *f"[  | 

*J (1S) 


From  the  shifting  theorem  already  presented, 

Y(f)  =  f %ft)  1  Jr 

=  Uro  Cjm 
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And  finally. 


Y(-f)  ~  Jf  £ j  c-6)  # Me t)J -  Q  (f)  H  ( -f)  (1? ) 

It  can  be  shown  similarly  that, 

¥  HCO]  *  f  C*)Mct)  (18) 

Clearly,  convolution  in  one  domain  Is  simple  multiplication 
In  the  other  domain.  Although  not  needed  for  the  pending 
development  of  the  discrete  Fourier  transform,  another 
Important  relationship,  known  as  the  correlation  theorem, 
can  be  appropriately  dealt  with  here.  The  correlation 
integral , 

/*  00 

2  (*)  »  j  Jto  f  O**.)  M*-  (19) 

has  an  operation  with  which  It  forms  a  Fourier  pair  as  did 
the  convolution-multiplication  operations.  This  theorem  can 
be  established  as  before, 

Z(£)  “J  (*•■**.)  *  Mt 

e»s  (3.4hf +)*!*•  t-j  f 

-  (20) 
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The  term  in  brackets  is  the  complex  conjugate  of  H(f)  and 
is  denoted  by  E*(f)  in  the  final  form  of  the  theorem  below. 


(»^  (t*  -t)  al'I'J  —  M 

(21) 

We  will  now  continue  on  to  the  discrete  Fourier 
transform  starting  with  a  continuous  waveform  h(t).  The 
waveform  is  sampled  or  multiplied  by  a  string  of  delta 
functions ,  s( t) . 

m 

t)s(t)  ~  £  (22) 


where  capital  delta  is  the  sampling  interval.  The  infinite 
sum  is  not  realized  and  must  be  windowed,  in  this  example, 
by  v(  t ) — 1  for  Mt  S(M-1)4  and  zero  elsewhere.  So  that 

now , 

M- 1 

^  Ci)  *  (*)  ^  T  (*-/*»*  )  (23) 

Recalling  the  convolution  theorem,  the  multiplications  in 
the  time  domain  correspond  to  convolutions  in  the  frequency 
domain  with  the  following  results.  H(f)  is  convolved  with 
the  window  functions'  spectrum  and  will  have  the  apparent 
effect  of  introducing  ripples  because  of  the  window's 
significant  sidelobes.  The  rippling  may  be  minimized  by 
choosing  a  window  function  with  small  sidelobes  at  the  cost 
of  other,  perhaps  more  acceptable,  spectral  degradation. 


E(f)  Is  convolved  with  the  sampling  function  and 
S( f )=I ( f-n/4  )  has  made  the  spectrum  periodic  with  respect 
to  the  Interval  7=1/4  .  The  spectrum  coming  from  a  real 
waveform  is  first  symmetric  about  f=0  as  shown  before.  Now 
because  of  Its  periodocity  the  spectrum  is  symmetric  also 
about  f=F/2  (the  Nyquist  or  folding  frequency).  One  final 
step  remains.  The  Fourier  spectrum  is  also  taken  at 
discrete  points  1/T  apart.  The  result  in  the  time  domain  is 
the  convolution  of  the  sampled,  windowed  signal  with 
I(t-nT),  which  is  a  periodic  signal  with  T  as  its  period. 

The  FFT  algorithm  used  in  the  preprocessors  developed  in 
this  thesis  uses  a  Cooly-Tukey,  base  two  algorithm  [15]. 

t 

This  is  documented  where  it  is  used  in  programs  included  in 
the  Appendices.  The  algorithm  uses  N  samples  where  N  is  two 
to  an  integer  power.  In  the  transformed  domain,  due  to  the 
symmetry  shown  above,  there  are  N  coefficients  (only  N/2  of 
which  are  unique  as  shown  in  figure  4).  The  original 
waveform  must  be  band  limited  prior  to  sampling  to  minimize 
aliasing.  The  resulting  frequency  spectrum  should  approach 
zero  at  the  folding  frequency  where  up  to  half  of  the  power 
can  be  aliasing  noise. 

In  the  following  two  chapters,  several  different 
preprocessing  algorithms  are  developed  and  their 
performance  compared  with  canonic  inputs.  In  Chapter  IV  the 


probleir  of  ship  identification  with  range  only  radar  is 
discussed  briefly  and  one  preprocessor  is  applied  to 
several  ship  profiles.  In  the  fifth  and  final  chapter, 
conclusions  are  drawn  from  the  work  done  here  and  follow-on 
efforts  are  recommended. 
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An  FFT  Spectrum 
Figure  4 
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II.  DIGITAL  l*ILL  IN  TRANSFORM 
BY  EXPONENTIAL  MAPPING 

In  the  first  chapter,  the  modulus  of  the  Pellin 
transform  was  shown  to  he  invariant  to  scaling.  A  detailed 
examination  of  the  mechanics  involved  suggests  an 
implementation  that  is  widely  used.  The  Mellin  transform  of 
a  t-domain  function  h(t)  is  given  in  equations  (7)  and 
O'). 

H<S)  *  £  J.i<) 

Hcs)  *  f  erV^  (s#) 

-4» 

Delta  has  been  added,  corresponding  to  the  sample  interval 
in  the  t-domain.  Again  it  is  noted  that  O')  is  a  Fourier 
transform,  where  s=-J2irm.  Solving  the  integral  for  the 
effect  of  a  t-scaling  by  the  factor  k, 

^  t'-'ott  *  ^  A  (*.  * 

=  «  sJ^A  His)  (24) 

Clearly,  in  the  Fourier  integral,  the  scaling  factor  k  has 
become  a  shift  for  which  the  modulus  of  the  transform  is 
invariant.  The  exponential  warp  alone  has  transformed  the 
scale  factor  into  a  shift.  A  prerequisite  is  that  the 
t-domain  signal  has  no  shift.  If  there  were  a  shift,  it 
does  not  transform  to  a  simple  factor  or  shift  in  the 
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Mellin  domain  and  so  cannot  subsequently  be  removed  by 
taking  the  magnitude  of  a  Fourier  transform. 

Implementation  of  a  discrete  Mellin  transform  is  as 
difficult  as  it  is  with  the  Laplace  transform  [13].  Once 
transformed,  characteristics  in  the  new  spectrum  are 
difficult  to  relate  to  the  original  signal  characteristics. 
One  hypothesis  relating  the  two  domains  is  generated  by 
comparison  to  the  Fourier  transform.  The  power  spectrum 
associated  with  the  Fourier  transform  can  be  used  tc  detect 
periodicities  in  the  physical  function,  since  the  wave 
numbers  at  which  sharp  peaks  of  the  spectrum  occur  give  the 
wavelength  of  such  periodicities.  3y  analogy,  the  positions 
of  the  peaks  in  the  spectrum  associated  with  the  Mellin 
transform  are  said  to  give  the  magnification  or  compression 
which  will  produce  features  in  the  physical  domain. 
Further,  this  stretching  and  compressing  is  identified  as 
periodic  in  nature  [21].  This  seems  unlikely  because  the 
Mellin  is  invariant  to  magnification/compression  and  does 
not  behave  well  (a  scale  factor  that  is  a  function  of  t  and 
k(t)  as  seen  in  (24).  The  Fourier  spectrum  models  the 
original  signal  by  a  set  of  weighted  sinusoids  of  varying 
frequencies,  and  therefore,  naturally  display  periodicity 
and  is  invariant  to  shift.  The  Mellin  is  also  a  weighted 
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sun  of  sinusoids,  but  whose  magnitudes  are  inversely 
weighted  by  t. 

rm 

flU>  ~  j  Act j  - — -  (25) 

Values  for  h{t)  for  0  t  1  are  far  more  important  to  the 
sum  than  those  beyond  that  point.  This  characteristic 
contributes  to  the  difficulty  of  realizing  discrete  Laplace 
and  Mellin  transforms  and  is  a  major  topic  covered  in  this 
chapter . 

The  numerical  approximation  of  the  Fourier-Pellin  (FI*) 
transformation  by  exponential  warping  Is  functionally 
diagrammed  in  figure  5.  A  FORTRAN  program  using  the 
algorithm  described  in  this  chapter  is  included  in  Appendix 
A.  Referring  to  figure  5,  the  input  samples  are  from  a 
pulse  whose  duration  is  finite  and  less  than  that  of  the 
sample  window.  For  this  case,  no  spectral  distortion  is 
experienced  by  filling  zeros  behind  the  sampled  data.  The 
only  effect  of  the  zero  filling  is  to  add  spectral 
resolution  in  the  frequency  domain.  Once  through  the  first 
FFT  block,  and  the  magnitude  taken,  Ef  is  symmetric  about 
zero  and  at  the  folding  frequency.  Thus,  for  N  time  samples 
and  filled  zeros,  there  will  be  N/2  unique  spectral 
coefficients.  This  unique  spectrum  is  interpolated, 
resampled  as  a  warped  signal  and  fed  to  the  final  FFT 
block.  A  correction  is  added,  and  the  modulus  is  taken.  The 
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Fourier-Mellln  by  Exponential  Warping 
Figure  5 


resulting  FM  features  are  invariant  to  shifting  and  scaling 
In  the  time  domain.  The  FIT  block  was  covered  In  sufficient 
detail  In  the  preceding  chapter.  Some  effort  will  be  spent 
in  discussing  the  warping  itself  and  the  need  for,  and  the 
development  of,  a  zero  point  correction. 

A.  ALGORITHM  DEVELOPMENT 

This  section  uses  its  own  notation  to  address  the 
requisite  exponential  sampling.  The  series  to  be 
transformed  is  h( f ) .  Its  Mellin  transform  is  E(m)  where  m 
is  the  Mellin  frequency  in  s=-j2frir.  The  transform  is, 

Hc~>  -  (26) 

Letting  f=Fexp[x]  as  before,  where  F  is  added  corresponding 
to  the  sample  interval 

H(~y)  =  J*  A  (re.*)  (2?) 

The  need  is  to  evaluate  M  equally  spaced  samples  in  x,  at 

Ot  X,  XX,  .  .  .  ,  <  x  (28  ) 

while  the  data  consists  of  N  equally  spaced  samples  in  f, 

0,F,zr,  ,  C  F  (29) 

Assuming  that  F  is  a  small  enough  interval  to  properly 
characterize  h(f),  it  is  sometimes  advisable  [9,10]  to 
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choose  the  exponential  sampling  interval  (X)  such  that  the 
largest  intersample  spacing  in  h(Fexp[x])  is  equal  to  kF 
where  k=l .  The  other  set  of  conditions  used  to  uniquely 
specify  the  new  samples  are  that  f=F  and  x=0  will  he  the 
lowest  limit  for  Interpolation,  while  (N-l)F  and  (P-l)X  are 
equated  as  the  upper  limit.  The  first  requirement 
constrains  the  choice  of  1“  by 


(30) 


Meeting  the  second  condition,  the  end  points  in  each 
sampled  series  are  equated  yielding, 

(H-ti  *  eCM")x  (31) 


Substituting  (31)  into  (30)  while  applying  the  exponential 
series  approximation  gives, 


X- 


jL 

H  -  I 


(32) 


One  more  substitution,  (32)  back  into  (31)  sets  up  the 
desired  result  where  M  is  now  specified  to  exponentially 
sample  from  f=F  to  f=(N-l)F  with  kF  being  the  largest 
interval  between  samples. 


M 


JL  ch-i) 

Jr 


(33) 
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As  N  gets  larger,  M  explodes.  If  N=l€  then  f*=41 ,  If  N=32 
then  M*106,  If  N=64  then  M=261,  and  so  on.  This  strict 
requirement  can  he  compromised  depending  upon  the 
application.  The  other  extreme  [18]  is  the  criterion  that 
requires  the  smallest  interval  between  samples  to  equal  the 
interval  between  uniformly  spaced  Nyquist  sampling,  with 
the  intervals  increasing  exponentially  thereafter.  The 
specification  permits  analyzing  frequencies  approaching  the 
largest  which  can  be  analyzed  with  uniform  spacing.  This 
greatly  reduces  the  number  of  samples  and  decreases  the 
required  processing  time,  but  is  limited  in  application. 
Two  factors  mitigate  the  stringent  requirements  Imposed  by 
(33)  where  k=l.  First,  using  N  unique,  uniform  samples 
yields  N/2  unique,  uniform  samples  in  the  FFT  domain.  A 
related  consideration  is  that  the  values  of  the  spectral 
components  approach  zero  at  the  folding  frequency  because 
the  original  signal  has  been  band  limited  and  some  over 
sampling  is  normally  recommended.  Secondly,  the  inverse  f 
weighting  apparent  in  equation  (25) 

i.  •  JLC-C)  .i#**^*^ 

j  e*  olf  (25') 

attaches  a  decreased  importance  to  h(f)  near  the  folding 
frequency  fn.  These  two  effects  combine  to  permit  a  much 
lower  sampling  rate  once  the  modulus  of  the  FFT  has  been 
taken.  In  spite  of  this,  the  stiffer  rule  (33)  is  used  for 
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this  work  to  generate  the  best  FM  domain  possible.  Whatever 
rule  is  used,  once  the  exponential  sample  points  have  oeen 
computed  for  an  FM  preprocessor,  they  needn't  be 
recomputed,  but  can  be  stored  for  rapid  access  during  the 
interpolation.  Some  interpolation  must  be  performed  to 
approximate  the  spectral  values  of  the  new  sample  points. 
Third  or  even  forth  order  Lagrange  polynomials  have  been 
recommended  and  used  for  this  purpose  with  apparent 
success  [9,10,19].  The  advantage  of  the  Lagrange  technique 
is  that  its  notation  is  particularly  compact,  consisting  of 
simple  summations  and  repeated  products  that  lend 
themselves  to  digital  realization.  Unfortunately,  for  the 
data  sets  used  in  this  thesis,  the  third  order  algorithm 
was  observed  to  behave  poorly,  adding  a  ripple  in  regions 
of  rapidly  changing  slope.  That  this  might  be  the  case  was 
suggested  by  the  advice  that  Lagrange  is  very  good  near  the 
central  data  point  when  the  order  of  the  polynomial  is 
known  to  be  the  same  as  the  order  of  the  approximation, 
otherwise  it  is  test  left  alone  [20,2lJ.  In  its  place,  a 
second  order  polynomial  was  used  to  interpolate  the  warped 
samples  with  results  that  were  nearly  indistinguishable 
from  the  actual  waveform,  as  seen  later  in  figure  8. 
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9.  ZEFO  POINT  CORRECTION 

Another  problem  becomes  apparent  when  the  Mellin 
transform  of  h(f)  Is  recalled. 


H(S )  » 

V 

&C4) 

''•40 

(34) 

-J  2 m . 

The 

exponentially  sampled 

wavefo  rm 

described  above  Is  applied  to  an  FFT  block.  As  f  approaches 
the  folding  frequency,  h ( f )  tends  to  zero.  Unf ortuna tely , 
as  f  approaches  zero,  the  value  of  h(f)  is  not  zero.  In 
fact  it  is  frequently  rather  high.  To  make  matters  even 
more  interesting,  the  left  integral  in  equation  (34^ 
clearly  shows  that  the  closer  to  zero  f  gets,  the  more 
important  h(f)  becomes  to  the  integral.  Several  solutions 
to  this  problem  are  considered  below. 

One  practical,  simple  approach  is  to  set  the  DC  (ie, 
f*0)  term  of  the  FFT  to  zero.  The  effect  is  nothing  more 
than  removing  a  DC  level  back  in  the  signal  domain,  but 
Mellln  transformation  into  the  FM  domain  leaves  the 
spectrum  dependent  upon  the  scale  factor  k  [22].  Setting 
the  f=0  coefficient  to  zero  corresponds  to  setting  h(f)  to- 
zero  for  0  f<l,  where  the  unity  upper  limit  is  chosen 


3? 


without  loss  of  generality.  The  resulting  transform  of  a 
scaled  signal  domain  h(f/k)  Is 

[  i  £*\ (351 

J 

which  is  obviously  dependent  upon  k.  The  closer  k  is  to 
unity,  the  smaller  the  effect.  By  increasing  the  f  spectral 
resolution,  the  error  can  be  reduced.  The  error  may  be 
insignificant  for  many  applications  [5-8],  but  the 
technique  should  be  used  with  care. 


The  first  solution  has  highlighted  the  need  for  a  zero 
point  correction.  Another  common  solution  [10-11]  is 
developed  by  breaking  the  integral  up  as  before.  Again 
using  unity  as  the  upper  limit  of  the  left  integral  while 
remaining  general  in  application, 

£ xcot'"jf  ♦ 

i  (3«) 


Two  assumptions  are  made  to  get  the  correction  term.  First, 
that  h(f)  remains  a  constant  h(0)  over  the  interval  of  f 
between  zero  and  one.  Second,  and  not  as  easily  accepted, 


..  V  u>  JL+  o 


=  e*®  =  o 


(37) 


Equation  (37)  pretty  well  shows  why  this  assumption  is 
suspect,  but  playing  along  for  the  momment ,  the  question  is 
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reserved  for  a  later  detailed  look.  Accepting  the 
assumption,  the  correction  factor  becomes, 

J?  fo) 


Z(s)  - 


.y  _  £<*L 


(38) 


Because  B(s)  Is  a  complex  function,  Z(s)  must  be  applied 
(added  to  the  Imaginary  part  of  the  succeeding  Fourier 
transform)  before  the  magnitude  is  take  to  remove  scaling 
dependence.  This  correction  Is  specifically  derived  for  use 
with  a  continuous  Fourier  transform  such  as  the  optical 
Fourier  In  Imaging  systems  with  the  added  stipulation  that 
h(f)  be  nearly  constant  over  the  range  0<f<k  where  k  is  the 
largest  scale  factor  expected.  If  an  FFT  is  employed  to 
make  this  final  transform,  another  correction  should  be 
applied  as  shown  by  Zwicke  and  Kiss  [11]  below.  This 
correction  factor  differs  from  the  first  due  to  an 
invariant  property  of  the  FFT.  The  FFT  of  two  unit  step 
functions  that  vary  only  in  scale  are  balanced. 


J*->*  £•*»' 


~  V  XtJt*  //vf 


(39) 


where  m  and  p  are  arbitrary  integers  greater  the  zero  and 
less  that  M.  Successive  FFT.  coefficients  are  summed,  and 
the  average  value  of  the  contributing  terms  taken  resulting 
in. 


r  -  J—  -  cc-6  ( it  J;  /m) 


(40) 
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This  is  then  multiplied  hy  h(0)  to  arrive  at  the  FFT  Mellin 
correction  factor. 


Z,„  CA)  * 


(41) 


When  k/M  is  small  the  imaginary  term  dominates  and  the 
correction  factor  approaches  that  used  in  the  continuous 
case  (38).  Most  of  the  work  done  for  the  thesis  on  the 
exponential  algorithm  was  done  using  the  inappropriate  zero 
point  correction  (38).  Since  its  discovery  was  coincident 
with  that  of  more  powerful  methods  discussed  in  Chapter 
III,  little  data  was  taken  using  (41). 


To  hound  the  error  involed,  an  acceptable  h(t)  is 
defined,  windowed  and  transformed  using  the  Mellin 
integral.  The  window  limit  is  then  allowed  to  grow 
unbounded  and  the  resulting  expressions  are  interpreted  as 
the  error.  It's  been  assumed  that  for  0<t<l,  h(t)=h(0). 
Since  the  error  resulting  from  the  assumption  in  equation 
(37)  arises  from  this  interval  alone,  t>l  is  ignored  for 
the  time  being.  Warping  the  signal  as  before  h(x)=h(0)  for 
all  x<0  and  h(x)=0  for  all  t>0.  Evaluating  the  integral  to 
a  finite  window  width  (T)  , 
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The  term  In  brackets  Is  the  magnitude  of  the  contribution 
from  the  region  of  integration,  -T<x<0.  Note  that  it 
involves  a  sin()/()  term.  The  effect  is  to  add  a  peak  of 
(T)h(0)  at  the  origin.  The  size  of  the  peak  depends 
directly  on  the  window  border  T.  Letting  T  approach 
infinity  raises  the  spike  at  ui  =0,  with  lesser  peaks  at 
^  =(2i+l )  ‘H’/T,  where  i  is  an  integer.  Fach  of  the  subpeaks 
has  a  magnitude  of  2h(0)T/(  (21+1)10  .  Substituting  into 
the  second  relation  yields  the  envelope  (in  brackets)  and 


phase  as  T  tends  to  infinity. 
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The  error  bound  in  brackets,  does  not  depend  on  the  sample 
rate,  or  the  size  of  the  window.  Any  approximation 
approaching  zero  will  have  the  same  bound.  Although  the 
magnitude  of  the  error  is  in  a  convenient  form,  the  phase 
is  indeterminate.  For  the  correction  to  be  applied,  the 
complex  addition  must  occur  prior  to  the  modulus  being 
taken.  This  cannot  be  done,  leaving  the  error  uncorrected, 
but  is  bounded  by  2h(0)/uJ  for  continuous  and  aperiodic 
discrete  Fourier  transforms.  For  the  FFT,  equation  (43) 
does  not  bound  the  error.  The  sum  of  1/i  does  not  converge 
as  i  tends  to  infinity.  At  any  point  on  the  FFT  this  sum  is 
present  due  to  the  apparent  folding.  The  error  itself  is 
not  unbounded  because  phase  differences  in  the  sum  of  the 


errors  at  any  point  may  result  in  the  envelopes  adding 
destructively,  reducing  the  actual  error.  Adding  the 
envelopes  is  an  unrealistic,  worse  case  approach.  The  FFT 
correction  (41),  can  not  he  compared  to  (43)  for  these 
reasons.  The  other  error  correction,  using  the  assumption 
in  equation  (37),  can  he  compared.  The  first,  setting  h(0) 
to  zero,  or  Just  plain  ignoring  the  3<t<l  interval  have  the 
error  function  hound  in  equation  (43).  Although  only 
differring  hy  a  factor  of  two  in  magnitude,  the  constant 
phase  is  arbitrary,  and  equivalent  to  setting  T=0J  that  is, 
assuming  h(0)=0  over  Z<t<l.  This  was  the  very  problem  the 
correction  was  developed  to  remedy,  hut  is  without  effect. 
Equation  (41),  the  correction  for  the  FFT  is  not  completely 
accepted  hy  this  author,  albeit  no  real  empirical  evidence 
have  served  to  verify  or  dispute  the  claim.  Suspicions  are 
raised  on  two  aspects.  The  indeterminate  phase  of  the  error 
(43)  in  the  continuous  case  arises  naturally  from 
approaching  the  t=0  limit  with  the  Mellln  integral.  This 
quality  is  conspicuous  hy  its  absence  in  the  FFT  error 
correction.  The  FFT  error  correction  is  computed  hy  summing 
the  FFT  coefficients  in  the  complex  plane.  The  average 
position  of  the  resulting  polynomial  is  the  error 
correction  term.  However,  the  FFT  coefficients  of  a  finite 
duration  signal  are  the  values  of  the  signal,  evaluated  at 
M  evenly  spaced  points  about  the  unit  circle  [3,23].  If 
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the  sequence  is  a  constant,  the  average  value  Is  the  origin 
of  the  complex  plane. 

The  zero  point  corrections  for  the  Mellin  transform  are 
unbounded  at  =0.  More  time  could  have  been  spent 
determining  the  best  applied  correction  for  the  specific 
case  at  hand,  but  direct  methods  are  developed  in  the  next 
chapter  that  obviate  the  need  to  employ  the  correction  at 
all. 

C.  TESTS  AND  RESULTS 

It  is  worth  admitting  at  this  point  that  the  results 
using  the  exponentially  warped  algorithm  to  achieve  a 
discrete  Mellin  transform  have  not  been  good.  More  recently 
developed  techniques  in  Chapter  III  greatly  surpass  the 
results  reported  in  this  section.  Although  much  of  the 
theory  used  to  improve  performance  In  the  following 
chapters  could  have  been  used  here,  this  was  a  preliminary 
attempt  that  was  later  abandoned.  The  FM  processor 
described  in  the  previous  section  was  built  using  FORTRAN. 
Appendix  A  is  the  documented  program.  This  section  will 
review  the  processing  with  actual  plots  of  the  signatures 
at  different  stages,  discuss  the  required  tests,  introduce 
the  testing  approach,  and  finally  intepret  the  results.  The 
functional  block  diagram  of  the  preprocessor,  figure  5,  is 
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near  the  beginning  of  this  chapter  and  could  be  profitably 
reviewed.  Figures  6  through  9  represent  the  step  by  step 
view  of  the  signal  processing,  where  the  signal,  a  test 
shape,  is  shown  in  figure  €.  The  waveform,  appears  as  an 
envelope,  and  is  drawn  with  vertical  lines  Indicating  the 
sampled  series.  Figure  7  is  a  picture  of  the  FFT,  in  its 
sampled  version  showing  its  symmetries.  Figure  8  is  a  very 
close  approximation  to  the  continuous  periodic  transform 
achieved  by  filling  zeros  to  obtain  the  requisite  spectral 
resolution.  A  "+”  on  the  transform  plot  indicates  an 
exponential  sample  point  interpolated  from  the  sixteen 
unique  points  in  figure  7.  The  warped  samples  are  sent 
through  the  FFT  block  once  more  with  the  result  shown  as 
figure  9.  Heavy  spectral  coloring  by  the  h(0)/ui  correction 
factor  is  evident.  Only  the  first  half  of  the  spectrum, 
from  zero  to  the  folding  frequency,  is  valid. 

Complete  scale  invariance  was  never  realized  although 
its  effect  was  greatly  reduced.  All  the  testing  done  on  the 
exponential  algorithm  was  an  attempt  at  achieving  and 
verifying  shift  and  scale  invariance.  Much  effort  was  spent 
structuring  the  tests  to  avoid  the  effects  of  processing 
noise  so  the  actual  algorithmic  characteristics  could  be 
determined.  Along  the  way,  requirements  arose  to  select  a 
suitable  interpolation  polynomial,  an  optimal  zero  point 
correction  and  other  system  improvements.  In  all  cases,  the 
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test  procedure  was  to  first  select  canonic  test  shapes. 
Squares  and  triangles  were  most  frequently  used,  being 
shifted,  scaled  and  combined  to  determine  preprocessor 
characteristics.  Scaling  was  most  frequently  by  a  factor  of 
two  or  less.  This  corresponds  to  an  aspect  angle  change  of 
60  degrees  from  the  unsealed  case.  For  many  tests,  care  was 
specifically  taken  to  exactly  scale  a  sampled  series 
instead  of  the  waveform  envelope.  The  variation  in  input 
signal  when  this  isn't  done  is  evident  when  considering 
figure  10.  For  the  envelope  shown,  the  sampled  series 
cannot  reconstruct  the  same  signal,  and  may  result  in 
feature  space  variation.  Two  feature  qualities  were 
monitored  in  each  test  to  determine  preprocessor 
performance;  insensitivity  to  shifting  and  scaling,  and  the 
ability  to  differentiate  between  canonic  classes.  These 
qualities  were  measured  by  visual  comparison  of  the  JM 
features,  by  computing  the  correlation  coef f lecients  and 
the  mean  squared  error  between  the  feature  sets  of 
differently  scaled  similar  test  shapes,  and  by  computing 
the  distribution  of  the  error  over  the  feature  space.  The 
latter  test  was  an  attempt  to  locate  feature  regions  of 
class  commonality,  and  regions  of  distinction  between 
canonic  classes.  This  approach  allowed  macroscopic  and 
microscopic  examination  under  changes  of  scale,  shift  and 
shape.  The  clustering  and  separation  qualities  are  data 
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dependent.  Since  no  ship  radar  video  data  was  used,  all 
observations  were  with  respect  to  canonic  classes.  No 
strong  groupings  were  detected. 


As  alluded  to  earlier,  the  results  for  the  exponential 
algorithm  were  less  than  satisfactory.  Test  shapes  were 
carefully  designed  to  minimize  sampling  effects,  and  to 
ensure  low  side  lobes  in  the  frequency  domain.  Pany  tests 
were  conducted  using  each  of  the  zero  point  correction 
methods  with  varied  shapes,  sample  rates,  and  spectral 
resolutions.  Classification  on  the  basis  of  signal  shape 
was  very  poor.  The  strongest  correlation  was  between  shapes 
of  common  duration.  Table  1  shows  a  typical  result  of 
comparing  a  rectangular  shape  and  a  raised  ramp.  The 
scaling  in  each  case  was  by  2  (60  degrees).  Equation  (35), 
Z=h(0)/ui  was  used,  but  the  others  offerred  little 
improvement.  Consistently,  the  strongest  similarity  was 
shown  between  shapes  that  had  the  same  sample  length,  vice 
shape . 


TABLE  1 


Canonic  Shape  Fourier  -  Mellin 
Feature  Comparisons 


a.  Peak  Correlation  Values 
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Ill .  BISECT  WELL  IN  TRANSFORMS 


The  last  chapter  developed  a  rrethod  of  obtaining  the 
Mellln  transform  by  exponentially  warping  the  signal  prior 
tc  using  an  FFT  block.  This  technique  Is  referred  to  as  the 
fast  Mellln  transform  (FMT).  Although  the  promise  of  scale 
Invariant  features  Is  attractive,  some  of  the  problems 
encountered  that  make  the  FMT  unattractive  are  reviewed 
here.  The  required  sample  rate  varies  with  respect  to  the 
data,  making  general  applications  difficult.  The  tendency 
Is  to  use  more  samples  than  required,  which  quickly  becomes 
costly  in  an  exponential  sampling  scheme.  The  need  to 
exponentially  warp  (to  interpolate)  a  set  of  new  samples, 
is  expensive  in  time  required.  For  true  scale  Invariance,  a 
correction  factor  is  required  but  because  of  the  Integral's 


unbounded 

natu  re 

at  zero  this 

correction  factor 

is 

indeterminant. 

Several  correction 

methods  have 

been 

employed, 

but 

they  depend 

on 

unspecified 

da  ta 

characteristics.  These  effects  combine  to  make  the  actual 
performance  of  the  algorithm  poor.  Although  scaling  effects 
are  mitigated,  they  remain  an  artifact,  which  is  disturbing 
to  classification  attempts.  This  chapter  outlines  the 
effort  to  remove  these  limitations.  Some  useful  Mellln 
properties  are  developed,  and  then  applied  to  establish 
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several  Direct  Mellin  Transforms  (EMTs)  which  were  built, 
and  their  performance  compared. 

A.  SOME  USEFUL  PROPERTIES 

Some  general  observations  are  made  here  about  the  Mellin 
transforms,  and  are  followed  by  some  specific  relationships 
which  were  derived  and  applied.  A  property  of  the  Mellin 
s-domain  is  that  it  is  unaffected  by  scaling  changes  in  the 
original  x-domaln.  Figure  11a  is  a  test  shape  in  the 
x-domain.  Two  features  are  identified  according  to  their 
amplitudes  A  and  B,  at  x=a  and  x=b  respectively.  The  ratio 
of  a/b  equals  c.  To  be  simply  scaled  by  k,  h(x/k)  must 
remain  the  same  in  all  aspects  except  that  the  distance 
between  features  has  been  changed  according  to  the  scale 
factor  k.  Figure  lib  shows  the  scaled  domain.  Features  A 
and  B  are  again  identified  at  their  scaled  positions  ka  and 
kb.  The  signal  property  maintained  by  scaling  is  relative 
positional  integrity,  that  is,  ka/kb  equals  c  as  before. 
The  positional  integrity  of  the  features,  the  ratio  of 
their  distance  from  the  origin,  to  that  of  another's  is 
unchanged.  Restated,  an  operation  0[h(x)]  in  the  x-domain, 
will  leave  the  s-domain  modulus  invariant  to  simple  scaling 
by  k  if 

0[A  <.*/*)]-  f<*M)  (44) 
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Note  that  the  entire  domain  is  scaled,  so  no  unsealed  shift 
in  the  domain  can  he  permitted  as  already  discussed.  In  the 
Beilin  s-domain,  any  simple  scaling  in  the  x-dorrain  results 
in  a  phase  distortion  (the  modulus  is  invariant  to  k). 
Manipulations  in  the  s-domain  will  leave  that  domain 
invariant  to  k,  as  long  as  the  modulus  is  modified  by  a 
multiplicative  factor  of  constant  phase.  That  is,  if  !&(s)! 
is  an  arbitrary  function  of  s  (except  that  it  does  not 
depend  on  k)  and  |M[h(x/k)]  !*'H(s) !  is  also  invariant  to  k, 
then  their  product  is  Invariant  to  k  and  the  x-domain 
remains  simply  scaled.  For  instance,  in  the  x-domain  the 
operator  0[h(x)]=x(h(x) ) ,  does  not  meet  condition  (44). 

7*  (45) 

So  the  Mellin  transform's  modulus  of  xh(x/k)  cannot  be 
invariant  to  k. 

In  Chapter  II,  a  error  was  apparent  because  h(x)  was  not 
equal  to  zero  for  x=0.  If  h(x)  could  be  modified  with  an 
operation  that  met  the  condition  of  equation  (44)  and 


always  produced  a  series  that  was  zero  at  x=0  a  general 
approach  can  be  developed.  Consider  two  operators. 


(46) 
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(47) 


Equation  (46)  will  produce  an  acceptable  f(x)  as  df/dx=0  at 
x*0.  Equation  (47)  will  always  produce  the  required 
condition.  To  see  the  frequency  domain  equivalents,  we  must 
assume  that  the  Mellln  integral  exists.  Integrating  by 
parts , 

|  |  *  |  Has>  I 

O 


I  j*  4nc  |  »  s  JmjtC«/4)«*~Ur  I 


~  I~*  **~U*  j  *  /*  (48) 

The  result  In  the  frequency  domain  is  consistant  with  the 
conditions  stated  above.  The  limit  of  h(x)  as  x  goes  to 
zero  or  to  Infinity  must  be  zero  if  equation  (48)  Is  to  be 
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true.  Similarly,  under  the  same  conditions,  it  can  be  shown 
that 


l*t[  £  (*  &*>*»]{  =  I  <**'>  “<J>I  (49) 

Equations  (48)  and  (49)  clearly  do  not  apply  where  h(0)  is 
not  equal  to  zero.  However,  by  using  the  x(d( h( x) )/dx ) 
modifying  operation  of  (47),  a  series  can  always  be  zero  at 
i*0.  The  function  h(x)  is  further  constrained  by  the  fact 
that  it  must  fall  off  faster  than  1/x.  This  assumption  must 
be  valid  and  the  modifying  operator  must  be  applied  in  the 
x-domain  for  the  Mellin  integral  to  exist  in  general.  A 
Mellin  transform  of  a  function,  after  having  a  modifier 
applied,  will  be  called  a  modified  Mellin  transform. 

f®  JA<«)  s 

“  Jo  (50) 

The  integral  is  close  to  a  form  which  is  realizable,  except 
for  the  upper  limit.  Tor  a  finite  sampled  series,  h(n)  will 
be  assumed  zero  outside  of  the  interval  0in<N.  This 
truncation  effects  the  transform  of  an  otherwise  infinite 
series.  In  this  application  the  Mellin  is  applied  to  an  FIT 
frequency  spectrum.  The  truncation  in  the  frequency  domain 
is  due  to  band  limiting  the  signal  prior  to  sampling. 
Scaling  in  the  time  domain  will  not  result  in  simple 
scaling,  but  will  add  a  dependence  on  the  scale  factor  Jc 
that  can  not  be  removed  by  the  transform.  An  approach  by 
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Prost  and  Goutte  is  used  to  predict  the  si2e  of  the  error 
[24,25].  First  a  suitable  function  will  be  selected  and  a 
relative  error  of  truncation  (RET)  determined  and  applied 
to  two  scalings.  The  relative  difference  of  the  feature 
space  is  found  and  identified  as  the  error.  Remembering 
that  this  is  applied  to  a  frequency  spectrum,  dh(f)/df  is 
approximated  by  a  function  of  the  form. 

J-lco  r  -f 

If  f  (51) 

The  modified  Beilin  transform  of  (51)  over  a  finite 
frequency  range  would  be  approximately, 


The  lower  limit  has  been  set  in  a  manner  to  be  consistant 
with  Plancherel's  theorem.  A  convenient  worse  case 
assumption  is  that  the  lower  limit  is  essentially  zero,  but 
this  depends,  in  general,  on  F  and  the  data  itself. 
Ha*(s,f)  converges  toward  Ea(s),  equation  (50),  in  the  mean 
square  .as  F  tends  to  infinity.  The  mean  square  error  will 
not  provide  an  expression  for  the  error  that  can  be  used  to 
correct  for  it  but  is  useful  to  measure  the  effect  of  the 
truncation  in  more  general  terms.  A  relative  error  of 
truncation  (RET)  can  be  computed  if  the  error  is  assumed  to 
be  distributed  evenly  over  the  range  from  which  Ha*(s,f)  is 
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computed.  By  employing  Plancherel's  theorem,  the  mean 
squared  values  are, 


<C‘ 

Jo  Jo 

s  CF  f*-  oi  f 

E  J0 

The  RET  Is  defined  and  solved  as, 

?6T  •  (~~T-  )'=  *  l)  * 

But  the  limit  F  depends  not  only  on  the  pass  hand  F,  but  on 
the  scaling  In  the  f-domaln.  A  relative  error  (RE)  between 
a  truncated  spectrum  and  a  scaled  and  truncated  signal 
would  be  more  complex,  but  worth  the  effort. 

+  r-n)-(  i/A)  *'-kF(zA'Lfr'+4F+i)  \ 
i  -  (zrztr-ti)  1 

where  0<k<l.  Two  observations  should  be  made  here.  First, 
the  relative  error  of  truncation  (55)  and  the  relative 
error  or  difference  between  two  truncated  spectrums  scaled 
differently  (56)  both  depend  on  F.  F  depends  on  the  cut  off 
frequency  of  the  low  pass  filter  prior  to  any  sampling.  It 
is  often  chosen  to  minimize  aliasing  depending  upon  the 


(54) 


(55) 
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limitations  of  the  sampling  circuit.  Second,  RI  depends  on 
k  as  well.  The  importance  of  scaling  differences  to  this 
data  type  can  he  readily  seen.  If  equations  (55)  and  (56) 
are  valid,  and  if  the  range  of  k  can  he  hounded  (a  design 
specification)  F  can  be  chosen  to  realize  a  stated  RE.  Or 
if  F  is  fixed  and  k  hounded,  the  RE  may  he  determined  for 
evaluation.  If  the  data  type  is  not  appropriate,  a  more 
representative  function  may  he  determined  and  used  in  place 
of  equation  (51)  to  attain  better  expressions  for  RET  and 
RI. 

3.  ALGORITHM  DEVELOPMENT 

Part  A  above  provided  some  background  for  making  the 
Direct  Mellin  Transforms  ( DMTs )  in  this  section.  In  this 
presentation  the  Justification  is  given  with  the 
application  first.  However,  chronologically  the  neat 
package  presented  above  was  preceded  by  extensive 
evaluation  of  empirical  results,  not  vise  versa.  Appendix  3 
documents  the  FORTRAN  implementation  of  all  the  algorithms 
developed  in  parts  1  and. 2  below.  Figure  12  is  a  functional 
block  diagram  of  a  generalized  TM  preprocessor  using  a  DMT 
to  show  the  simplicity  with  which  it  can  be  applied,  as 
opposed  to  the  FMT  covered  in  the  second  chapter,  figure  5. 
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1 .  First  Difference  Approximations 

Although  developed  through  a  different  rationale, 
this  first  algorithm  was  developed  by  Zwicke  and  Kiss  [ll]  . 
Starting  with  a  sampled  series  hi,  the  series  is  operated 
on  by  the  modifier  defined  In  equation  (51),  using  the 
first  backward  difference  to  approximate  the  derivative 
with  respect  to  x.  Unit  step  size  is  assumed. 


/* 
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(57) 


Taking  the  trapezoidal  rule  to  evaluate  the  modified  Mellin 
integral  (52)  while  recalling  that  h(0)=0,  and  h(N)  is 
assumed  zero. 


where  s*-J2fhm/M.  The  complex  coefficients  are 


SH*-  C0S  2 IT  (* /m  )  -ji  si*  Zft  0*0 


They  can  be  calculated  off  line  and  stored  to  produce  Just 
the  desired  characteristics.  The  factor  (2-^m/M)  could  be 
any  number  that  produces  an  Interesting  feature.  Undesired 
features  need  not  be  computed  (what  in  general  was  a  N  by  H 
matrix*  where  M  Is  the  number  of  Fellin  transform 
coefficients  and  N  the  number  of  x  sample  points).  If  a 
relatively  small  number  of  features  Is  required,  perhaps 
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the  processing  will  he  manageable.  Notice  that  no  zero 
point  correction  is  required.  Only  data  changes  contribute 
to  the  transform.  These  observations  are  valid  for  any  of 
the  modified  Direct  Wellin  Transforms  developed  in  this 
section, 


By  using  a  central  difference  instead  of  the 
backward  difference,  a  simillar  result  is  obtained. 


IV'l 

Hfct  <V»v)  *  Zj  ^'*♦1  ^  /z 


(60) 


Other  numerical  integrations  may  be  used  with  improved 
results,  and  other  methods  can  be  used  to  increase  the 
order  of  the  approximation. 


To  test  the  algorithms,  a  ramp  and  inverse  ramp  were 
used  in  a  scaled  and  unsealed  mode.  The  ramps  and  their 
scalings  are  shown  in  figure  13.  Figure  14  is  the  analytic 
results  of  both  waveforms  plotted  with  the  transform  found 
using  equation  (60).  This  is  a  dramatic  improvement  over 
anything  used  with  the  methods  discussed  in  the  previous 
chapter.  The  noise  of  the  signal  appears  to  be  diverging  as 
frequency  grows,  but  over  the  range  plotted,  the  appearance 


is  that  of  an  algorithm  trying  to  do  well. 


2 .  Second  Difference  Approximations 

A  second  difference  algorithm  can  be  achieved  by 
proceeding  as  before.  The  pure  scaling  operator  used  to 
prepare  the  signal  is, 

OtA<*) ]z  ^  -ZJL  +  zA  ) 


And  the  new  algorithm  is, 


N-l 
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s  *  I 
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(62) 


Other  second  order  operators  have  been  used,  and 
partitioned  into  forward  and  backward  difference  variations 
to  (61),  but  this  appears  to  be  a  basic  and  useful  form. 
The  color  l/(s+l)  is  present  to  approximate  the  modified 
Mellin  of  equation  (60).  The  term  (s+1)  is  valid  assuming 
that  dh(x)/dx  is  exclusively  upper  bounded  by  ln(x)/x  as  it 
approaches  zero  or  infinity.  For  comparison,  another  second 
difference  algorithm  was  developed  based  on  the  modifier 
(x(d/dx) (x(d/dx)h(x) . 


M-f 


(63) 


This  is  roughly  the  sum  of  the  methods  defined  in  equations 
(60)  and  (62)  above.  The  assumption  for  deriving  the  color 
1/s  is  even  less  restrictive  than  before,  but  the  algorithm 
performs  poorly  compared  to  (62).  These  transforms  depend 


6? 


on  second  difference  characteristics,  but  are  modified  by  a 
1/s  term  which  has  a  stabilizing  effect.  This  is  equivalent 
to  a  division  by  x  and  integration  in  the  x-domain.  The 
order  of  the  approximation  has  been  increased  by  the 
modifier.  Results  using  equation  (60)  and  (62)  should  be 
alike,  figure  15  is  the  results  of  using  (62)  compared  to 
the  closed  form  solution  to  the  figure  13  test  shape 
transforms.  An  improved  performance  over  (60)  is  seen  over 
some  of  the  range,  but  a  drop  as  frequency  increases 
degrades  the  accuracy  in  figure  15b. 


3.  Higher  Ilfference  Approximations 
Higher  difference  approximations 
For  instance,  one  algorithm  depending 
difference  is. 


can  be  developed, 
upon  the  third 


(64) 


The  performance  of  higher  order  algorithms  becomes 
increasingly  suspect  because  of  the  extreme  weighting  they 
apply  to  different  parts  of  the  series.  Different 
algorithms  exist,  but  this  weighting  is  always  a  factor.  A 
smoother  transform  is  achieved,  but  a  large  error  is  likely 
to  develop  due  to  the  algorithm's  dependence  on  higher 
order  derivatives  and  the  nature  of  sampled  data.  However, 
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these  comments  are  speculative  since  they  were  not 

confirmed  experimentally. 

4.  Higher  Order  Integrations 

Higher  order  integration  rules  should  he  able  to  be 
used  with  a  corresponding  improvement  in  performance.  One 
using  the  first  difference  with  Simpson's  rule  was 
implemented  with  di ssappointing  results.  Figure  16  is  the 
result  of  such  an  implementation.  The  droop  for  the  ramp 
input  is  apparent  even  though  not  present  in  the 
trapezoidal  rule  used  in  subsection  1  above.  The  higher 
frequency  error  is  also  more  prevalent  than  before.  A 
program  error  is  of  course  suspected,  but  was  never  found. 
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I-MELLIN  F 

Simpson's  Rule 
Figure 
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IV.  CLASSIFICATION  PR2PR0CISS I NG 


The  problem  of  determining  important  signal 
characterises  can  be  approached  in  at  least  two  different 
ways.  First,  by  trying  to  learn  what  is  important  in  human 
recognition,  and  then  trying  to  adapt  a  machine  to  emulate 
that  behavior.  Or  second,  by  using  successive  transforms  to 


remove  information 

known 

to  be 

superfluous 

to 

classification,  while 

keeping 

eno  ugh 

inf orma tion 

to 

reliably  assign  an  object  to  a  class.  Addressing  the 
former,  even  though  it  is  difficult  to  determine  specific 
details,  some  key  aspects  of  human  visual  recognition  are 
dlscernable.  Chief  among  these  is  that  the  intensity  level 
of  a  scene,  or  object,  does  not  appear  to  be  as  important 
as  the  relative  position  of  the  edges,  i.e.  the  shape 
separating  different  intensities  and  frequencies  [25-281. 
Examples  in  scene  analysis  show  clearly  that  the  edges  or 
shapes  are  far  more  critical  to  human  recognition  than  the 
relative  power  differences  themselves.  The  invariant  shapes 
or  angles  in  scene  analysis  find  their  analog  in  ratios 
between  similiar  points  in  differently  scaled  time  series. 
Examination  of  many  range  only  radar  video  ship  signatures 
has  provided  a  basis  for  noting  that  the  relative  position 
of  time  domain  features  remain  constant  over  changes  in 
aspect  angle,  while  the  relative  intensity  of  the  features 


vary  greatly  as  shown  in  figure  17.  As  the  aspect  changes, 
if  this  set  of  ratios  remains  constant  within  a  ship  class, 
then  the  set  is  identified  as  information  worth  preserving 
for  the  classifier.  Conversely,  since  relative  intensltiy 
is  not  a  stable  measure  over  aspect  angle,  or  from  among 
different  ships  of  one  class,  that  information  should  be 
intentionally  removed  to  provide  tighter  natural 
clustering,  with  the  minimum  number  of  features. 

A.  INFORMATION  REQUIRE!  TO  CLASSIFY 

There  are  two  preconditions  that  must  both  exist  for  a 
set  of  possible  input  signatures  to  separate  into  distinct 
classes.  First,  the  features  of  a  particular  class  must 
have  some  common  characteristic  about  them,  and  second, 
this  characteristic  must  in  some  way  be  unique  with  respect 
to  other  classes.  The  assumption  in  existing  radar 
signature  classification  projects  is  that  there  is  enough 
information  in  the  signatures  to  permit  this 
classification.  Short  of  actually  trying  to  classify  with  a 
set  of  realistic  signatures,  the  analytical  determination 
that  sufficient  information  is  present  in  a  set  of  all 
possible  ship  signatures  is  difficult  to  approach. 

To  establish  how  well  the  autonomous  classifier  is 
performing  some  measure  of  the  classifiabili ty  of  the  set 
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received  signals  is  desired.  Failing  this,  some  discussion 
of  the  Information  capacity  of  the  preprocessor  should  at 
least  be  considered.  The  preprocessor  produces  the  FFT 
magnitude  of  a  sampled  signal  as  the  output  of  the  first 
stage.  Most  of  the  unique  positional  relationships  of  the 
signal  upon  which  human  recognition  apparently  depend  has 
been  destroyed.  Next,  the  Mellin  transfrom  stage 
effectively  distorts  the  signal  and  uses  the  magnitude  of  a 
second  Fourier  transform  as  the  output  features.  Signals 
reconstructed  on  the  basis  of  FFT  phase  information  alone 
usually  provide  sufficient  similarity  to  be  associated  with 
the  original  signal,  whereas  reconstruction  on  the  basis  of 
magnitude  does  not  retain  any  significant  detail  except 
when  the  signal  is  symmetric  [29].  The  data  is  also  masked. 
For  most  applications,  the  data  is  frequency  shifted, 
filtered  and  sampled  as  a  baseband  signal.  This  windowing 
masks  the  magnitude  characteristic  and  is  specifically 
designed  into  the  processing.  The  resulting  FM  features  are 
insensitive  to  positional  and  scaling  relationships  that 
are  necessary  in  human  visual  recognition. 

Analytic  support  is  also  available  to  quantify  the 
importance  of  relative  position  of  events  [29].  By 
considering  rms  error  (due  to  spectral  phase  and  amplitude 
quantization  for  random  signals)  it  has  been  concluded  that 
approximately  two  more  bits  are  required  for  the 
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quantization  of  phase  information  than  amplitude  for  the 
same  rms  error.  A  separate  analysis  applied  distortion  rate 
theory  to  real-part,  imaginary-  part,  and  magnitude-phase 
encoding  of  the  DFT  of  random  sequences.  The  result  was 
that  phase  required  1.4  hits  more  storage  than  magnitude 
for  a  similar  error  [30].  A  third  approach  concluded  that 
the  Fourier  phase  includes  1.6  hits  more  information  than 
the  magnitude  [31] .  This  was  based  on  analysis  of  image 
reconstruct  ion  from  kinoforms  (phase-only  holograms).  The 
fact  that  phase-only  reconstruction  preserves  much  of  the 
correlation  between  signals  would  suggest  that  the  location 
of  events  tends  to  be  preserved.  Further,  it  seems  that 
this  information  is  lost  by  taking  only  the  magnitude  of 

the  Fourier  transform.  Another  interesting,  albeit 

informal,  view  is  apparent  as  one  considers  the  phase-only 
signal  as  a  spectral  whitening  process. 

W  =  /F(f;|e->9'  (e5 

For  reconstruction  by  phase  alone,  where  the  magnitude  is 
set  to  one! 

[Fciii  (66) 

Since  the  received  radar  signature  will  have  an  abundance 
of  low  frequency  spectral  lines  and  smaller  high  frequency 
components,  the  low  frequency  information  is  not  weighted 
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as  heavily  as  the  higher  frequency  information  in  the  phase 
reconstructed  signal.  This  seems  like  it  would  accentuate 
sharp  changes  in  the  reconstructed  signature  without 
removing  relative  location  information.  The  result  is  the 
summation  of  the  different  frequency  components,  all  with 
zero  phase  (i.e.,  no  positional  or  amplitude  distribution 
information  can  remain). 

Considering  the  Mellin  transform  next,  in  a  continuous 
case,  the  exponential  warp  does  not  lose  any  information. 
To  zero  the  first  data  sample  as  required  by  the  Chapter 
III  modifications,  will  surely  destroy  information,  but 
this  may  be  confined  to  the  EC  term  alone.  The  information 
lost  during  the  final  transform  and  magnitude  is  difficult 
to  assess.  Increased  masking  occurs  due  to  the  spectral 
truncation,  so  actual  information  loss  may  not  be  as  great, 
but  masking  distortion  may  be  greater  than  before.  So 
approximately  two  bits  of  information  are  lost.  Only  a 
quarter  of  what  was,  remains.  Information  is  also  lost  when 
the  transforms  are  normalized  in  the  processing  so  that  any 
power  calculation  is  also  meaningless.  Some  interesting 
questions  arise.  After  removing  positional  relationships, 
scaling,  and  power,  what  signal  qualities  remain  and  are 
they  useful  in  classification?  Although  it  is  true  that 
this  insensitivity  may  add  a  certain  robustness  to  the 
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system,  the  arbitrary  loss  of  valid  classification 
information  should  be  minimized. 

By  examining  the  quality  that  must  be  ignored,  and  by 
comparing  its  removal  to  what  is  actually  removed  by  the 
processing,  an  interesting  result  will  develop.  The  effect 
of  a  time  domain  shift  on  the  frequency  domain  is  an 
additive  phase  term,  linearly  related  to  the  frequncy  of 
the  coefficient,  as  in  equation  (2).  most  of  the  structure 
of  the  signal  is  held  in  the  phase  relationships  with 
respect  to  the  fundamental  and  higher  frequency  terms.  So 
shift  can  be  defined  as  the  phase  of  the  fundamental 
complex  coefficient.  By  setting  the  phase  to  zero,  and 
adjusting  the  other  coefficients  according  to  their 
component  frequency,  the  structure  of  the  signal  is  not 
lost  but  reconstructed  about  the  fundamental  as  before.  If 
there  are  N/2  spectral  phase  angles,  only  the  fundamental 
needs  to  be  zeroed  to  remove  the  shift.  If  the  information 
is  contained  uniformly  in  the  spectral  phase,  then  only  2/N 
of  this  structural  information  needs  to  be  removed.  When 
the  magnitude  is  taken  to  produce  shift  invarient  features, 
all  of  the  phase  relationships  are  destroyed.  (N-2)/N  of 
the  Information  once  held  in  the  phase  was  removed 
needlessly.  The  amount  of  information  lost  removing  the 
shift  can  be  made  arbitrarily  small.  As  N  grows  unbounded, 
the  amount  of  information  that  needs  to  be  removed  tends  to 
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zero.  This  result  seems  to  be  supported  by  experience.  With 
enough  samples  of  a  shape,  its  position  with  respect  to  the 
observation  field  is  immaterial.  In  theory  the  principle  is 
sound,  but  some  practical  limitations  may  degrade  predicted 
performance.  Recalling  that  an  exponential  warp  translates 
scaling  to  shifting  generalizes  the  result  a  bit  further. 

Ad/Jt)  =  JL  (6?) 

The  same  priciple  that  permits  simple  shift  removal  is  also 
valid  for  the  removal  of  scaling  dependence  as  well.  Using 
the  Pellin  transform,  scaling  dependence  may  be  removed  by 
zeroing  the  fundamental  and  adjusting  all  the  other 
coefficients  as  described  above.  For  the  FM  preprocessor, 
the  information  lost  removing  the  scaling  and  shifting 
dependence  may  be  made  arbitrarily  small  by  increasing  the 
number  of  spectral  samples  used.  Since  the  number  of 
spectral  samples  can  be  increased  by  filling  zeros  onto  the 
finite  signal  in  the  original  domain,  this  process  does  not 
effect  the  data  sample  rate.  This  approach  was  not  verified 
experimentally,  but  represents  a  potentially  powerful  tool 
to  analyse  and  improve  the  extracted  feature  space. 


79 


B .  RANGE  CMY  RADAR 


This  section  addresses  ship  classification  on  the  basis 
of  information  gathered  from  a  range  only  radar  video  ship 
signature.  An  example  of  such  a  signature  has  already  been 
considered  as  figure  16.  Classification  by  range  only  radar 
signatures  is  subject  to  the  same  distortions  discussed 
above.  Typically,  the  radar  return  is  detected  and  isolated 
in  a  range  gate  that  is  sampled  and  digitized.  The 
rectangular  sampling  window  can  be  considered  the  range 
gate  itself.  The  range  gate  is  designed  to  ensure  that  the 
included  range  is  greater  than  the  maximum  ship  length  so 
that  the  time/range  windowing  has  no  effect  on  the 
frequency  spectrum  of  the  signature,  other  than  increased 
spectral  resolution.  The  placement  of  the  ship  signature  in 
the  window  is  not  set,  neither  from  encounter  to  encounter, 
nor  from  pulse  to  pulse  (Jitter).  The  total  effect  joins 
together  to  produce  the  framing  distortions.  Signature 
scaling  results  from  viewing  the  ship  from  different  aspect 
angles.  The  sampling  rate  must  be  done  at  more  that  twice 
the  inverse  of  the  resolution  of  the  receiver.  Quantization 
levels  are  chosen  in  a  manner  to  reduce  that  predictable 
random  noise  to  an  acceptable  level.  The  pulse  to  pulse 
jitter  is  an  ever  present  characteristic  of  the  radar 
problem,  but  at  a  normal  resolution  (greater  than  26  feet) 
integrating  the  return  in  the  time  domain  removes  the 
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Jitter  effect,  reduces  scintillation,  and  Improves  the 
resolution  of  the  signature.  Although  predetection 
(coherent)  integration  is  more  efficient,  post  detection 
(noncoherent)  integration  is  more  common  because  of  the 
convenience  of  not  having  to  preserve  the  radar  frequency 
(FF)  phase.  For  post  detection  integration  of  n  pulses,  the 
signal  to  noise  ratio  would  be  something  less  than  n  times 
the  signal  to  noise  ratio  for  one  pulse  [32].  More 
important  to  the  recognition  problem  itself,  for  a  stable 
system  by  the  law  of  large  numbers  [33],  fluctuation  of  the 
average  value  of  the  return  will  be  overcome.  That  is,  with 
the  integration  of  n  pulses  the  resolution  (R)  will  become 
finer  as 


"R  a  ij  5  (68) 

with  a  probability  of  1-E ,  where  S  squared  is  the  variance 
of  the  signal  from  pulse  to  pulse.  For  very  high 
resolution,  the  cost  of  making  the  signature  stable  with 
respect  to  the  integrator  becomes  prohibitive.  It  has 
become  convenient  to  integrate  the  spectrum  of  the 
signature  because  the  jitter  effects  can  be  completely 
removed.  Another  limiting  factor  in  integrating  over  a 
period  of  time  is  that  the  position  and  aspect  cf  the 
target  are  dynamic.  They  change  with  time.  A  conceptually 
attractive  solution  is  a  recursive  filter  which  weights  the 
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Integrated  pulses  such  that  the  older  they  become,  the  less 
weight  is  accorded  to  them. 

If  the  course  and  speed  of  the  target  ship  are  known 
from  measurement  of  the  target  track,  it  is  possible  to 
infer  the  aspect  of  the  ship.  The  range  profile  can  give 
some  estimation  of  size.  Unfortunately,  the  three 
dimensional  change  in  aspect  angle,  commonly  suffered  by  a 
ship  presents  more  than  just  a  video  signature  scaling 
change.  The  radar  cross  section  of  even  individual 
structural  components  of  the  radar  target  changes  with 
respect  to  the  aspect  angle.  The  composite  effect  is  that 
ship  signatures  vary  greatly  with  aspect  angle.  The  radar 
is  an  electromagnetic  sensor,  reacting  to  energy  reflected 
from  the  target.  These  reflections  are  a  result  of 
scatterers  that  are  related  in  dimension  to  the  wavelength 
cf  the  illuminating  energy.  Because  of  the  great  difference 
in  wavelength  between  light  and  microwaves,  what  can  be 
"seen"  by  radar  may  be  duite  different  than  that  seen  by  an 
eye.  Also,  when  measuring  size  or  any  distances  with  a 
radar  of  high  resolution  (less  than  £0  feet),  an  error  can 
be  incurred  since  the  extremities  of  the  target  are  not 
always  good  scatterers  .  Echoes  from  the  forward  or  stern 
portions  of  the  target  might  be  observed  in  the  noise, 
especially  for  a  relatively  low  power  radar  [32].  After 
reviewing  hundreds  of  signatures,  it  appears  that  major 


82 


features,  such  as  overall  ship  length  and  dominant  mast 
structure  are  frequently  discernable,  but  vary  in  relative 
amplitude.  Resonance,  shadowing  (one  reflector  hiding 
another),  multipath  returns,  the  mapping  of  three 
dimensional  aspect  changes  onto  a  one  dimensional  time 
series,  and  the  amplitude  and  phase  of  component  returns 
summing  constructively  or  destructively  to  cause  a 
scintillation  of  the  composite  target.  Some  of  the 
variations  caused  by  these  conditions  can  be  lessened  by 
Integration  but  major  effects  remain  causing  the  signature 
to  vary  ir.  shape  and  content  with  the  aspect  angle.  For 
this  reason,  a  class  feature  volume  cannot  be  reduced  to  a 
single  point,  but  will  remain  a  hypervolume  in  the  feature 
space  even  In  the  ideal  case.  Any  selection  of  features 
should  try  to  minimize  this  volume.  Features  should  be 
selected  that  are  relatively  insensitive  to  known 
superfluous  effects. 

C.  CLASS  DISCRIMINATION 

In  the  last  chapter,  the  major  concern  was  removing  two 
sources  of  variance  with  no  classifying  value.  Algorithms 
were  developed  and  canonic  shapes  generated  to  verify  the 
algorithms  and  demonstrate  the  invariance  to  scaling.  In 
the  same  manner,  this  chapter  has  reviewed  information  loss 
and  process  masking.  The  question  of  whether  sufficient 
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information  is  present  to  classify  was  raised.  To  answer 
this  question  some  simple  classes  of  canonic  figures  are 
defined,  and  put  through  the  entire  FM  preprocessor 
documented  in  Appendix  B.  The  LMT  algorithm  found  to  offer 
the  best  scale  factor  rejection  in  Chapter  III  was  used  to 
generate  the  final  features.  The  algorithm  chosen  is  based 
on  a  second  difference  modification  and  is  defined  in 
equation  (62), 

W-/ 

Csj  s  s  ^  (62) 

/n=i 

After  canonic  tests  were  made,  preprocessor  performance  on 
several  ship  signatures  was  recorded.  Although  this  was 
premature  in  the  logical  test  sequence,  the  results  are  of 
some  interest. 

1 .  Test  Shapes  and  Results 

Four  test  shapes  were  used.  Figure  18  shows  the  test 
shapes.  All  were  scaled  and  shifted  originally  to  test  for 
algorithm  verification  and  demonstrate  scale  invariance.  In 
this  series  of  tests  they  were  left  fixed  and  used  in 
different  combinations  to  try  to  detect  shape  presence  in 
the  F^  feature  space.  The  object  is  to  differentiate 
between  different  canonic  classes.  Figure  19  shows  a 
comparison  of  the  shapes,  rectangle  and  triangle.  A  test 
combining  the  rectangle  and  triangle  was  the  subject  of 
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figure  20.  Finally  in  figure  21  a  rectangle  with  two 
triangles  Is  shown  In  the  FM  feature  space.  It  is  clear 
froir  the  plots  that  most  of  the  scale  variance  has  been 
removed.  Just  as  important,  some  quality  does  remain  that 
differentiates  between  the  canonic  shapes.  A  square  in 
general  can  be  differentiated  from  a  triangle.  A  "ship" 
with  a  single  mass  of  superstructure  can  be  separated  from 
one  with  two  such  masses. 

Although  it's  clear  from  the  plots  that  there  is  a 
unique  quality  left  in  the  feature  space  to  allow  the  time 
domain  shapes  to  be  classified,  some  quantified  measure  of 
system  performance  is  required.  The  magnitude  of  the 
feeature  vectors  have  all  been  normalized  with  respect  to 
the  first  coefficient,  so  in  that  region  little 
discrimination  can  be  expected.  For  higher  Beilin 
frequencies,  noise  dominates.  A  region  of  coefficients, 
11-100  was  chosen  to  classify  the  shapes  by  correlation. 
The  results  are  included  as  Table  2.  The  improvement  in 
performance  over  that  shown  by  Table  1  in  Chapter  II  is 
dramatic.  The  methods  in  that  earlier  test  resulted  in  the 
observed  length  of  an  object  being  the  distinguishing 
criterion  for  classification.  Using  methods  supported  in 
Chapter  III,  variations  due  to  scaling  and  shifting  of  the 
original  domain  have  been  removed.  The  features  now  reflect 
the  shape  of  the  object  in  the  time  domain.  An  unusual 
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S-MELLIN  FREQUENCY 


Shape  1  FP  Features 
Figure  20 
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APE  2 


1 


S-MELLIN  FREQUENCY 


Shape  2  FM  Features 
Figure  21 
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effect  is  noted  with 
analysis  in  Table  2b. 
correlated  to  TRI1  than 
the  reverse.  The  results 
has  greatly  simplified 
canonic  shapes  above. 


respect  to  the  difference  squared 
Although  TRI2  is  more  closely 
RECT2  the  squared  error  shows  just 
are  encouraging.  The  preprocessor 
the  classification  problem  for  the 


2.  Ship  Signatures  and  Results 

A  single  ship  was  used  to  make  these  preliminary 
tests  for  this  thesis.  Signatures  were  taken  every  ten 
degrees  around  a  ship  from  zero  to  fifty  degrees.  The 
results  are  plotted  and  compared  over  twenty  degree  aspects 
in  figures  22-23.  The  signatures  are  the  result  of  very 
high  resolution  radar  signature  data  that  has  been  degraded 
and  smoothed  to  a  lower  resolution  with  essentially  no 
noise  present.  Recalling  that  the  purpose  of  the 
preprocessor  was  to  remove  variance  due  to  pure  shifting 
and  pure  scaling,  leaving  enough  information  for 
classification,  to  the  eye  there  seems  to  be  little 
encouragement  from  these  results.  It  is  recalled  that  a 
goal  of  this  preprocessor  is  to  make  the  classifiers  Job 
easier  by  removing  dependence  on  shifting  and  scaling  of 
the  original  data.  The  information  that  remains  depends  on 
unspecified  signal  characteristics  that  here  appear  to 
useful  in  discriminating  shape  classes  and  possibly  ship 
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classes.  However  no  real  conclusion  can  be  drawn  at  this 
point  because  of  the  small  data  base  and  the  absence  of  an 
automatic  classifier  to  generate  an  optimal  feature  space. 
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Table  2 


Canonic  Shape  Fourier  -  Mellin 
Feature  Comparisons 

a.  Peak  Correlation  Values 
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RECT/2 

TRI 

RECT 

1.00 

0.95 

0.50 

BECT/2 

- 

1.00 

0.52 

TRI 

- 

- 

1.00 

TRI/2 

- 

- 

- 

uared  difference 

between  features. 

RECT 

RECT/2 

TRI 

RECT 

.  000 

.011 

.010 

RECT/2 

- 

.  000 

.015 

TRI 

- 

- 

.  000 

TRI/2 

— 

— 

— 

TRI/2 
0.42 
0.41 
0.98 
1 .00 


TRI/2 
.011 
.014 
.011 
.  000 
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V. 


The  preprocessor  design  began  by  considering  a  generic 
classification  systeir.  A  distinction  was  drawn  between  the 
classifier  and  the  preprocessor.  The  preprocessor  is 
problem  specific.  It  assists  in  the  classification  by 
extracting  a  set  of  features  for  the  classifier.  The 
extracted  feature  space  has  enhanced  natural  clustering  on 
the  basis  of  shape  by  removing  information  that  was 
extraneous  for  classification.  Two  useless  characteristics 
were  identified  as  shifting  and  scaling.  The  preprocessor 
was  designed  to  remove  dependence  on  these  two 
characteristics  by  using  the  invariant  properties  of  a 
Fourier  transform  followed  in  series  by  a  Mellin  transform. 
The  resulting  set  of  coefficients  are  Fourier-Mellin  (FM) 
features . 

A.  REVIEW 

In  Chapter  II  a  Mellin  transform  was  developed  using  the 
conventional  digital  processing  approach  which 
exponentially  warps  the  domain  and  then  transforms  the 
spectrum  by  an  FFT .  This  method  is  sometimes  identified  as 
a  fast  Mellin  transform  (FMT).  The  unbounded  behavior  of 
the  exponentially  warped  frequency  spectrum  was  shown  to 
result  in  a  function  that  could  not  be  transformed.  That 
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is,  the  warped  function  has  no  transform  because  the  Mellin 
integral  is  indeterminate  at  the  lower  limit.  Error 
correcting  techniques  cannot  compensate  for  the  effect 
because  the  error  itself  cannot  be  computed  in  general.  The 
bound  for  the  error  was  found  and  seen  to  be  the  envelope 
for  existing  error  correction  functions. 

Chapter  III  developed  some  useful  properties  of  the 
Mellin  transform  that  were  used  to  modify  the  signal  so 
that  the  pitfalls  isolated  in  Chapter  II  were  avoided.  This 
was  done  by  modifying  the  input  to  the  Mellin  transform  to 
always  be  transformable.  To  simplify  the  implementation  and 
to  control  the  effects  of  sampling,  a  direct  Mellin 
transform  was  used  for  the  development  of  the  modifiers. 
Several  suitable  modifiers  were  determined  and  tested  with 
differently  scaled  inputs  for  which  the  closed  form 
solution  was  known.  The  direct  Mellin  algorithm  that 
produced  features  closest  to  the  closed  form  solution  was 
chosen  for  use  in  the  preprocessor.  With  the  modifications 
in  place,  the  preprocessor  was  tested  and  shown  to  produce 
features  that  were  invariant  to  shifting  and  scaling.  It 
was  also  shown  that  the  features  retained  enough 
information  to  classify  canonic  shapes. 

Chapter  IV  discussed  what  type  of  information  is 
required  for  classification.  Signal  structure  or  shape  was 
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identified  as  key  inf  orirat  ion .  A  discussion  of  what 
information  is  required  for  classification  and  a  means  of 
keeping  the  signal  structure  intact  throughout  the 
preprocessor  was  advanced ,  hut  not  empirically  verified. 


B.  FUTURE  WORK 


The  design  FM  preprocessor  does  produce  a  feature  space 
with  enhanced  clustering,  but  problems  remain  to  be 
resolved  before  the  full  potential  of  the  system  can  be 
realized.  There  are  three  extant  conditions  that  detract 
from  the  performance  of  the  implemented  preprocessor. 
First,  the  preprocessor  is  not  computationally  efficient. 
Second,  most  of  the  signal  structure  that  should  be  vital 
to  classification  is  obviously  lost.  And  third,  a  complete 
verification  of  performance  has  not  been  conducted.  The 
current  preprocessor  produces  an  enhanced  feature  base  for 
a  classifier,  but  attention  to  these  main  weak  points  will 
greatly  improve  the  applied  techniques. 


1.  Efficient  Processln 


A  direct  Mellin  transform  using  N  spectral  samples 
to  transform  into  M  Mellia  spectral  coefficients  requires  M 
by  N  multiplications.  If  only  a  few  coefficients  are  to  be 
used  then  the  number  of  multiplications  may  be  small.  Using 
the  FMT  requires  about  N(20+ln(N))  arithmetic  operations 
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for  the  interpolation  and  the  FFT.  If  twenty-five  or  more 
Mellin  coefficients  are  required,  the  FMT  is  faster.  The 
modifiers  used  to  prepare  the  spectrum  for  the  direct 
Mellin  transform  will  also  work  for  as  FMT,  but  this  should 
be  demonstrated  experimentally.  Because  the  FMT  directly 
weights  the  lower  numbered  samples  it  may  produce  more 
accurate  results  as  well. 

2.  Conservation  of  Information 

Chapter  IV  discussed  the  type  of  information 
required  for  visual  pattern  recognition  in  humans.  The 
preservation  of  this  information  should  be  a  specified 
design  goal  for  the  preprocessor.  The  preprocessors  built 
for  this  thesis  removed  much  of  these  vital  signal 
properies.  A  means  was  introduced  to  limit  or  control  the 
loss  of  structural  detail  by  zeroing  the  fundamental  phase 
and  adjusting  each  of  the  remaining  complex  coefficients  to 
reconstruct  phase  relationships.  This  may  be  done  in  the 
frequency  domain  as  described,  or  by  a  similar  operation  in 
the  time  domain,  shifting  the  centroid  to  zero.  In  either 
case,  to  transform  more  information  about  the  signal  will 
effectively  increase  the  sensitivity  of  the  features  to 
characteristics  in  the  time  domain.  This  increase  in 
sensitivity  needs  measured  to  confirm  the  approach.  It  is 
also  possible  that  continued  selective  zeroing  of 
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V 


coefficients  may  offer  improved  performance  or  robustness 
to  the  system  as  a  whole. 

3 .  Verification 

Although  the  improved  performance  was  demonstrated 
with  respect  to  ealier  Ftf  digital  preprocessing,  a  direct 
improrement  factor  needs  to  be  established.  Time  domain 
correlation  should  be  used  as  a  measure  of  original  signal 
classifiability.  This  measure  needs  to  be  compared  to  the 
FF  domain  correlation  recorded  as  Table  2  in  Chapter  IV. 
Next,  the  preprocessor  or  an  improved  version  will  have  to 
be  married  to  a  classifier  and  realistic  data  used  to 
evaluate  its  effect  on  the  classification  system.  The 
preprocessor  built  was  design  to  be  used  on-line.  An 
on-line  classifier  needs  to  be  built  as  well. 

The  systematic  evaluation  of  ship  profiles  using  FI* 
features  is  still  a  requirement.  For  several  ships,  FH 
features  must  be  singled  out  and  plotted  together  as  a 
function  of  aspect  angle.  The  purposes  are  to  establish  a 
range  of  aspect  over  which  classification  may  be  possible, 
to  evaluate  changes  of  structural  content  as  discussed  in 
section  two  above,  and  to  determine  the  beam  signature  as  a 
classification  "node".  Although  these  purposes  assist  in 
the  system  design,  the  third  is  more  of  an  operational 
necessity.  All  ship  signatures  will  degrade  to  the  beam 
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aspect  "node"  so  that  this  hypervolume  in  the  feature  space 
is  occupied  in  common.  Therefore  the  beam  condition  must  be 
detected  and  withheld  from  ever  entering  the  classifier. 
Nor  should  beam  signatures  be  used  for  classifier  training. 
The  beam  "node"  needs  to  be  determined  separate  from  the 
classifier.  Initially  this  information  can  be  plotted  to 
examine  feature  behavior  and  confirm  the  methodology, 
automated  methods  will  quickly  follow.  These  analysis 
functions  may  never  reside  in  the  classification  system 
itself,  but  must  be  a  part  of  the  tools  used  to  develop  a 
workable  classification  system. 


100 


nnnn  nnonnnnoonnonr>nnnnnnnr>nnnnnn 


APPENDIX  A 
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I************************************************ 

THIS  PBOGBAM  IS  A  DIGITAL  IMPLEMENTATION  OP 
A  FAST  FOOBIEB  TBANSFOBM  PQLLONED  B7  A  MELLIN 
T BANS FOBS .  THE  DIGITAL  IMPLEMENTATION  OF  THE 
MELLIN  IS  AS  AH  EXPONENTIALLY  SAMPLED  3PECTEUM 
THEN  BOH  TBBOOGH  AH  FPT  AND  COBBECTED  FOB  THE 
EBBOB  BY  ANY  ONE  OF  SBYEBAL  COBBECTIOH 


SOBBOOTINES,  COBCTX 

******************* 


****** ***************************** 


BECOBD ED  PLOTS  /  PLOT  HOMBEB  OSING  BECOBD  CALLS 
-TIME  FUNCTION  (CONT)  /  1 
-SAMPLED  TIHE  FUNCTION  /  1 
-FFT  (CONT)  /  2 

-EXPONENTIALLY  SAMPLED  FFT  /  2 
-UNIFOBHLY  SAMPLED  FFT  /  3 
- DISCRETE  MELLIN  FEATOHES  /  4 


BECOBDED  PLOTS  /  PLOT  HOMBEB  OSING  BECANG  CALLS 

-time  function  jconti  /  i 

-FFT  ONIFOBMLY  SAMPLED  /  2  (MAG  &  PHASE) 
-MELLIN  FEATOBES  /  3  (HAG  &  PHASE) 


150 

200 


DIMENSION  XBEAL  (200)  ,  XIMAG  (200)  ,XINT(200)  ,T(64)  , 

1 FBT  (200) 

I  SCALE  =  31 
BU*5 
H*2**H0 
NO*  7 
B*2**NO 

CALL  SAMP  (XBEAL* XIMAG. H) 

CALL  BECOBD  (XBEAL, XIMAG ,N,ISCALE) 

COMPOTE  THE  ACTUAL  SAMPLE  POINTS  AS  PBOVIDSD  BY  THE 
SAMPLED  YIDEO. 

NO»5 
N-2**NO 


CALL  SAMP  (XBEAL.XIMAG, N) 

CALL  BECOBD  (XBEAL,  XIMAG  ,  Hf  ISCALS) 

ISCALE  *  31 
10*7 
l»2**L0 
ro  200  1*1, L 
IFiX.GT.  Njj  GO  TO  150 


XiJt(I) 

(I)! 


PBT  ( 

GO  T 
XIHT 
PBTJI) 
CCHTIB 


0 

ill 


XIMAG 
200 
)  *0.0 
*0.  0 
OE 
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310 


300 

C 

C 

C 

C 

c 

c 

c 

c 

c 


400 

C 

C 

C 

c 

c 

c 

c 

c 

c 

399 


401 


415 

410 

C 

C 

c 

c 

c 

c 

c 

c 

ccc 


CALL  FFT (XINT,PRT,L,LU) 

CALL  FFT  (XBEAlt, XIMAG,  N  ,NU) 

CALL  BECOBO  (XaEAL.  XIMAG.  N, ISCALE) 

CALL  BECOBO  (XINT,PBT,L,  ISCALE) 

ISCALE  *  31 
00  300  I*  1 ,  N 

XBBA1  (X)  *SQRT  (XBEAL  (I)  **2+XIMAG  (I)**2) 

XXHAG  (I]*0. 

CONTINUE 

CALCULATE  T-HE  MSB  SAflPLE  POINTS  ABO  IBTEBPOLATE 
TO  FIBO  THE  EXPOBEBTIALLX  SAMPLED  7 ALOES. 

CALL  HUPTS  (XINT,M,N) 

AFTEH  N EH  XHTEBVALS  CALCULATED  ABO  STOBEO  TEMPORARILY 
IB  7ECT0B  IINT.  STOBE  FOB  LATEB  PBIBTXBG  LB  T  VECTOR. 

DO  400  1*1  ,M 
PBT  (I)  *0.0 
(I)*X 


flHOE 

HITH  THE  BEH  TIMES  IN  VECTOR  XINT,  AND  THE  COBBEBT 
SPEC! HUH  SAMPLES  IN  VECTCB  XBEAL.  COMPUTE  THE  NEH 
EXPONENTIAL  SAMPLES  AND  ENTEB  THESE  INTO  VECTOR 
XINT  BT  USING  THE  CHOSEN  INTERPOLATION  METHOD. 
FINALLY.  BECOBO  THE  FIRST  SAMPLE  TO  BE  USED  LATER 
TO  CORRECT  THE  MELLIN  TRANSFORM  FOB  LOH  FREQUENCY 
LOST  IN  THIS  EXPONENTIALLY  SAMPLED  TECHNIQUE. 

FO-XBEAL(I) 

CONTINUE 

CALL  INTP2  (XBEAL, N, XINT ,M) 

BRITEJ2  ,40 1}  M 
FORMAT  (14) 

CALL  SMAX(XINT,PBT,M, SCALE) 

00  410  I*1,M 
XREAL  (I)*XIBT  (I) 

XINT  (II  "SCALE  *  XINT  (I) 

XIMAG  (I)*0 

HBITEJ2.4 151T (I) .XINT (I) 

FORMAT ($10. 5, 7X, FI  0.5) 

CONTINUE 

SUBMIT  THESE  EXPONENTIALLY  SAMPLED  VALUES  TO  THE 
FINAL  FFT  BLOCK. 

CALL  FFT (XBEAL, XIMAG, M, MU) 

APPLY  THE  CORRECTION  TERM,  FIND  THE  MAGNITUDE 
OF  THE  BOH  SCALE  AND  TIME  INVASIENT  FEATURES. 

CALL  COBCT1  (XBEAL, XIMAG, H,F0) 

CALL  COBCT2  (XBEAL, XIMAG, M,F0J_ 

CALL  BECOBO  (XBEAL, XIMAG , M, ISCALE) 

STOP 

END 
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nnoo 


C 

c 

c 

c 

c 

c 


ft***************************************** 

*  SOBBOUTINE  PPT#  GIVEN  THE  COMPLEX  * 

-  (2**HU»B||g  HILL  BETOBN  THE  * 


SAMPLES 


* 

*  N  CQEFPIfiCIENTS  USING  A  DPT.  * 

****************************************** 


102 


101 


100 


103 


SUBBOOTINE  PFT (XBEAL, XIMAG, N ,N0) 

DIMENSION  XBEAL  (N)  ,  XIMAG  (N) 

N2-N/2 

N01*NO-1 

R»0 

DO  100  L*1,N0 
DO  101  I»l!N2 
f«IBIIB  (K/2**N01,NU) 
ABG»6.2B3185*P/PLOAT 
C*COS 4ABG) 


(») 


S»SIN jABG) 

K1*K+1 

K1H2=»K1+N2 

1BBAL3XBEAL-(K1N2)  *C+XIMAG  (K1N2)  *S 
TIMAG«XIMAG(K1 N2J *C -XBEAL  K  1 N2) *S 
- I1N2)*XBEAL(K1)-TBEAL 


XBEAL 

XIMAG  |K1N2|.>XIMAG^ 
XBEAL  (K 1)  *X8EAL(Kl 
XIMAG  (K 1 j  = XIMAG (K1 
K*K+1 
K3fC  +  N2 

IP (K.LT.N)  GO  TO  102 
R>d 

N0l*N0l-1 
N2*N2/2 
00  103  K«1.N 
I-IBITB  (K-1.NO)  M 
IPil.LE.K^GC  ' 


K1J-TIMAG 
♦TBEAL 
♦TIM AG 


TBEAL»XBE2 


TO 
-(K) 
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TIMAG«XIMAG  K 


XBEAL 
XIMAG 
XBSAL 
XIMAG  , 
CCNTINO 
BETOBN 
END 


)3XBEAL(I) 
J  *  XIMAG  (I) 
J-TBBAL 
‘TIM  AG 


K 

K 

h 


PONCTION  IEITH 
PONCTION  IBITB (J, NO) 

J 1*  J 
IBIIB«0 
DO  200  1*1,  NO 
J2-J1/2 

IBITB*IBITB*2+  (J1-2*J2) 
200  J1*J2 

BETOBN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


10 

25 


30 


60 


C 

C 

c 


c 

c 

c 


c 

c 

c 


******************************************************** 

*  INTP2  IS  A  SECOND  OBDEB  INTBBPOLATION  BASED  ON  A 

*  PCLINOEIAL  EQUATING  TO  PIBST  AND  SECOND  DEBIVATZ VES 

*  APPROXIMATED  BX  CENTBAL  DIFFEBEHCES.  THE  INPUT  VECTOB 

*  <X>  CONTAINS  UNIFOBMLX  SAMPLED  DATA  HITH  UNIT! 

*  BETNEEN  CONSECUTIVE  SAMPLES.  VECTOB  <T>  IS  INPUT 

*  HITH  THE  NEH  SAMPLE  TIMES  AND  OUTPUT  HITH 

*  THE  NEH  SAMPLES.  THEBE  ABE  N  ONIFOBH  SAMPLES  IN  <X> 

*  AND  M  HABPEC  SAMPLES  IN  <I>. 
******************************************************** 


SOBBODTINE 
DIMENSION  X 


I 


CHOSE  THE  X  SAMPLE  TIME  CLOSES  TO  THE  HABPED 
TIME  HELD  IN  <X> 

DO  60  1 1* 1 # M 
ET  *  100 
DO  10  I Xs 1 , N 
IXTIME  *  IX  -  1 
ETABS  »  ABSJDT) 

CIST  »  FLOAT  (IXTIME)  -  Y  (IT) 

DIST  »  ABS(DIST) 

IF  (DTABS  .LT.  DIST)  GO  TO  25 
DI  *  Y(II)  -  IXTIME 
CONTINUE 

IXTIME  -  IXTIME  -  1 

J1  »  -1 

DO  30  J*1 ,3 

J1  =  J  ♦  IXTIME  -  1 


cg^tJc 


X(J1) 


NUE 

TIME  »  I(IY) 

YJII)  a  rCN- (X3 1 TIME , IXTIME) 

CONTINUE 

BETUBN 

END 

FCH  IS  THE  INTEBPOLATION  BOLE. 

FUNCTION  FCN(X,T,I) 

DIMENSION  XJ3) 

TX  a  FLOAT  (I) 

COMPUTE  THE  COEFFICIENTS. 

A  a  (X  (3)  -  2.*  1(2)  ♦  Xfin  /  2. 

E  =  jx?3)  -  Xlin  /  2.0  -  2.*  A  *  TX 
C  *  X(2)  -  A  *  IX** 2.  -  B  *  TX 

COMPUTE  THE  INTEBPOLATEC  VALUE. 


FCH  «  i  *  (1**2)  ♦  B  *  I  ♦  C 

BETUBN 

END 
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************************************************* 

*  SOBBOOTIME  IMTEBP  OSES  A  LAGBANGE  THIBX) 

*  CBDEB  METHOD  OB  A  SECOMC  OBDEB  POLYNOMIAL 

*  TC  COMPOTE  THE  INPOT  SAMPLE  MA7EF0BM. 

* 

*  THE  SOBBOOTINE  INPOTS  ABE: 

*  X-SAMPLED  SA7EFOBM 

*  M-THE  NOMBEB  OF  X  SAMPLES 

*  I-INPOT  SAMPLE  TIMES 

*  -OOTPOT  IIEBPOLATED  SAMPLES 

*  M-THE  NOMBEB  OF  I  SAMPLES 
************************************************* 


SOBBOOTINE  INTEBP  (X  ,  N.Y  ,M) 

DIMENSION  X4  (4)  ,  X  (N )  ,  X  (  M)  ,S  HFTX  ( 1 56) 

FIBST  CHOOSE  XTIME  NEABEST  EACH  I TIME. 

THEM  COMPOTE  THE  IMTEBPCLATED  7 ALOE  AT  THAT  PT. 


CO  100  1*1, N 
I 5*1+5 

IFII.GE.  (H-5) )  GO  TO  101 
SHFTX(I5)  *X  (I) 

GO  TO  100 
15*1-  (N-5) 

SHFTX  (l5)  *X(I) 

COMTINOE 
CO  40  1*1, M 
COMTINOE 
CO  60  II* 1 ,  M 
YP*100 

CO  10  IX* 1 ,M 
IXT-IX-1 
YPABS*ABS  (IP) 

DIST*ABS  (IU-I  (II)  ) 

IF  (IP ABS. LI. BIST)  GO  TO  25 

IP*Y  (IT) -IXT 

IX4*IX+5 

COMTINOE 

114*1X4-2 

CO  30  J*1 -4 

X4(J)*SHPTX(IX4+J) 

COMTINOE 

IM1*IX4+1 

10*1X4+2 

IP1*IX4+3 

IP2*IX4+4 

I  (II)  *ILAGB  (X4,  IP) 

COMTINOE 

BETOBN 

END 


nnnnnnn  o  cn  nnnn 


1 


FONCTION  XL IN  HAKES  A  LINEAH  INTBBPOLATION . 


FONCTION  YLIN (X4, YP) 
CIHEJ1SI01I  JW(«J 
IF(YP.LT.O.)  60  TO  50 
YLIN*X4  (2)  ♦YP*(X4(3)-X4  (2)) 
B  ETOBN 

0  YLIN*X4(2)-YP*<X4(1)-X4{2)) 

BETOBN 
END 


EDUCTION  YLAGH  COMPOTES  THE  LAGBANGE  MOLTI-PLIEES 
AND  HAKES  THE  I NTEH EOLATION  FOB  A  CHOSEN  OFFSET 
f BOH  THE  CENTEAL  SAHPLE  X4(2) 


C 

C 

C 

C 

C 

C 

C 

C 

C 


100 


FUNCTION  YLAGB (X4,YE) 

DIMENSION  X4 (4) 

CH1»-YP*  ( YP-T)  *  (YP-2) /6 . 0 
CO*  <YP**2-1)*  (YP-2)  /2.0 

YLAGB*CH1*X4  (1)  ♦C0*X4  (2)  ♦CP1*X4  (3)+CP2*X4  (4) 

BETOBN 

END 


******************************************************** 

*  SOBBOOTINE  NDPTS  CALCOLATES  THE  H  EXPONENTIAL  SAHPLE  * 

*  EGINTS  FBOH  THE  N  ONIFOBH  SAMPLES  OF  THE  EXISTING  * 

*  SPECTBOH  IN  PBEPABATION  FOB  AN  INTEBPOLATION.  THIS  * 

*  EXPONENTIALLY  SAMPLED  SET  OF  POINTS  ABE  STOfiED  IN  * 

*  THE  INPOT  VECTOB  X.  * 

******************************************************** 


SOBBOOTINE  NUPTS (X, B, N) 
DIMENSION  X(156) 
QNaFLCAT (N) /2 . 0  ♦  1.0 
EH=FLCAT(H)-1.0 
CELZ-LOGONJ/EM 
DO  100  X«1.jf 
SI» FLOAT (11-1.0 
“  '  SI* DEL! 


VALOE 


fINOE 
BETOBN 
END 
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c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


10 


20 

50 

C 

c 

c 

c 


*  mis  IS  A  STOB  PROVIDING  k  TEST  SKIM  RETURN 

*  TC  SIMULATE  k  SKIN  RETURN.  NOBS ALL I  THIS  SUBROUTINE 

*  HILL  ACTUALLY  SAMPLE  A  SKIN  RETURN. 

********* *********************************************** 

SUBROUTINE  SAMP (XB, XI, N) 

DIMENSION  XB(N)  ,XI  (N) 

DESIGNATE  SCALE  FACTOR  AND  UNSCALED  TIME  SHIFT. 

SCALE  =  16.  /  32. 

SHIFT  *  0.0  /  32. 

SCALE  =  1.0/SCALE 
TO*  .5  ♦  SHIFT 
N1=N-1 

CO  100  1=1, N 
XB  ~ 

XI 


100 


1  VW  A* 

(I)  =0. 

(I)=0. 


CALCULATE  THE  TIME  OF  THE  SAMPLE. 

SK=  (FLOAT  (IJ-I.O)  /FLOAT  (N1) 

EUILC  THE  TEST  SKIN  RETURN. 

IN  THIS  (C  2-PER)  CASE  A  DOUBLE  PERI MID. 
T  SC  ALE  =  (SK-TO)  *SC ALE 
IHFM  =  1 

GO  TO  10 
GO  TO  50 


GO  TC  20 

.£)  *  10.0 


IF  (IHFM 
IF (IHFM 
H  =  8.  / 

IF  (TSCALE  .LT.  -H)  GO  TC  100 
IF  (TSCALE  .GT.  +  H)  GO  TO  100 
IF  JTSCALE  .LT.  TO)  - 

XB  (I)  =  (TO  -  TSCAtl 
GO  TO  100 

XB  (I)  »  H  *  10.0  -  (TSCALE  -  TO)  *  10.0 

GO  TO  100 

CONTINUE 

THE  FOLLOHING  HAVE  FORM  IS  A  SQUARE  HAVE  16  SAMPLES 
HIDE  CENTERED  AROUND  THE  SIXTEENTH  SAMPLE. 

SCALING  AND  SHIFTING  DONE  ABOVE  HILL  EFFECT  THE 
HAVEFCRM  ACCORDINGLY. 


TSCALE  =  (SK-TO) *SCALE 
EDGE  =  1.  /  4. 

IF  (TSCALE  .LT.  -EDGE) 

IF  JTSCALE  .GT.  +EDGE) 

XB (I)  =  1.0 

CONTINUE 

RETURN 

END 


GO  TO  100 
GO  TO  100 
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T| 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


******************************************************* 

*  THIS  SUBROUTINE  RECORDS  A  SELECTED  DATA  SET, 

*  AND  SCALES  THE  INDEX  TO  32  SAMPLES  (0-31). 

~  THE  INPOTS  ARE  XREAL.  XI MAG ,  AND  N  (THE 


.»  **««<*,  anu  n  iiaa 

ISC  ALE  IS  A  CONTROL 


50 


40 


60 


75 

100 


C 

C 


«  SOMBER  OP  SAMPLES) 

*  VARIABLE. 

*  IF  ISCALE  *  0  NO  SCALING  IS  PERFORMED 

*  SCALE  -  1,  AXIS  =  31 

*  SHIFT  =  1  (STARTS  AT  T  =  0  ) 

*  ISCALE  <  0  SHIFTING  AND  SCALING  OCCURS 

*  1 ISCALE!  =  AXIS  LEN .  A  SHIFT  IS  INCORP- 

*  SEATED  SO  THAT  AXIS  STARTS  AT  1. 

*  ISCALE  >  0  AXIS  SCALING  OCCURS  BUT  THERE  IS 

*  NO  SHIFT  (IE.,  THE  AXIS  STARTS  AT  0.) 

***************************** ************************** 

SOBROOTINE  RECORD  (XR.XI.N, ISCALE) 

DIMENSION  XB(N)  ,XI(N)  ,RlC  (200) 

SCALE*1 .0 
AXIS  a  31 
SHIFT  »  1 
H BITE  (2,50)  N 
FORMAT (14) 

IF  (ISCALE  .EQ.  0)  GO  TO  40 
AXIS  *  FLOAT (ISCALE) 

AXIS  *  ABS  (AXIS) 

IF  (  ISCALE  .LT.  0)  SHIFT  =  0 
CALL  SMAX  (XB, XI, N, SCALE) 

DO  100  1*1, N 

0) 


100 


SI=  (I-S HI f!)  *AXIS/( FLOAT  (Nl  -  1 .  ( 

** 2+XR  (I)  **2) 

CONTINUE 

HEC  (I)  =SCALE*R£C  (I) 

«  RITE  (2.75)  SI  .EEC  (I)  „ 

FORMAT (PI  0.5, 7X,P  T 0.5) 

CONTINUE 


RETURN 
END 


SUBROUTINE  SMAX  (XR,  XI,  N  , SCALE) 
DIMENSI CN  XB(N)  ,XI  (N)  ,T  (200) 
XMAX=1.0 
DO  100  I=T.N 

T  (I).=SQBT  'XR(I)  **2+XI  (I)  **2) 
IF(T(I)  .Li.XMAX)  GO  TO  100 


XHAX  =  T (I) 
SCALE=1 .0/1 
CONTINUE 
RETURN 
END 


/XMAX 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


90 

100 


c 

c 

c 

c 

c 

c 

c 

c 


90 

100 


************************************************** 

*  THIS  CORRECTION  SOBEOUTINE  OSES  ONE  Of  > 

*  THE  SIMPLER  C0BBEC1I0MS  FOB  THE  MELLIN  ' 

*  TBANSFOBM .  THE  COBBECTICN  IS  A  PURE  IMAGINARY. 

* 

*  COBBECTION  *  -FO/OMEGA 

*  THEM  A  MODIFICATION  IS  MADE  TO  THE  ENTIBE 

*  TRANSFORM  BX  A  MULTIPLICATION  BY  OMEGA. 

************************************************** 


SOBBOOTINE  C0BCT1 (IB, XI , N ,F0) 
DIMENSION  XB(N)  ,XI<N) 

DO  100  1*1,  N 

OMEGA  *  FLOAT (I)  -  1. 

XB  (I)  =  XB  (I)  *  OMEGA 

IF  (I  .EQ.  1).GO  TO  90 

XI  (I)  =  1X1  (I)  -  FO/OMEGA)  *  OMEGA 

GC  TO  100 

XI  (I)  =  FO 

CONTINUE 

BETOBN 

END 


******************************************************* 


*  THIS  CORRECTION  APPLIES  THE  MOBE  COMPLECATED 

*  EXPRESSION , 

*  CORRECTION  *  FO/2  ♦  JCOT  (FO/OMEGA) 


*  AND  THEN  MODIFIES  THE  ENTIBE  TRANSFORM  BY  1/OMEGA, 

*************************•**************************: 


*** 


SOBBOOTINE  CORCT2 (Xa,XI ,N,FQ) 

DIMENSION  XB  (N)  ,  XI  (  N) 

CO  100  1-1,  N 

OMEGA  =  FLOAT  (I)  -  1. 

XB  (I)  =  (XB  (I)  *  FO/2.)  *  OMEGA 

If  (I  .EQ.  1)  GO  TO  90 

XI  (I)  =  1X1(1)  -  (FO/2.)*COTAN  (OMEGA)  )*OMEGA 
GO  TO  100 

XI (I)  *  FO/2. 

CONTINUE 

BETOBN 

END 
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APPENDIX  B 


******************************************************** 

*  THIS  PROGRAM,  POURIER  DIRECT  MELLIN ,  TAKES  AN  INPUT 

*  HAVEFORM  FRCM  LOGICAL  DEVICE  2,  PERFORMS  AN  FFT 

*  FOLLOMED  BY  A  MDMT  CUTPUTING  THE  FEATURES  TO  LOGICAL 

*  DEVICE  3  FOR  LATER  PLOTTING  BY  MELPLT.  THE  MAJOR 

*  DATA  STRUCTURES  AND  SUBROUTINES  ARE  LISTED  BELOH 

*  IN  ORDER  OF  THEIR  APPEARANCE.  THE  SUBROUTINES  ARB 

*  DESCRIBED  IN  MORE  DETAIL  HHEBE  THEY  ARE  ACTUALLY 

*  LISTED  IN  THE  PROGRAM. 

******************************************************** 


HAJOB  DATA  STRUCTURES: 

< HFM>  -  THE  INPUT  HA7EFOBM  (REAL) 

<XFM>  -  THE  MELLIN  TRANSFORM  (COMPLEX) 

<CPHI>  -  REAL  MELLIN  COEFFICIENTS 
<SPHI>  -  IMAGINARY  MELLIN  COEFFICIENTS 
<STAND>  -  AN  ARRAY  HELD  FOR  LATER  COMPARISON 
OR  OTHER  USE. 

<PRT>  -  AN  ARRAY  NORMALLY  USED  TO  HOLD  REAL 
DATA  TEMPORILY.  A  WORK  SPACE. 

<IXT>  -  X  AXIS  TITLE  FOR  PLOTTING 
<IYT>  -  Y  AXIS  TITLE  FOR  PLOTTING 
<KEY>  -  NUMBER  OF  PLOT  THIS  GRAPH 


SUBROUTINES: 

HAVE  -  BEADS  AN  ARRAY  FROM  LOGICAL  DEVICE  2, 
AND  FILLS  ZEROS  TO  MAKE  A  TOTAL  OF  256 
SAMPLES.  THE  OUTPUT  IS  IN  <HFM>. 

FFT  -  AN  FFT  BLOCK 

COEF  -  COMPUTES  THE  MELLIN  TRANSFORM  SAMPLE 
HEIGHTS.  THESE  ARE  COMPLEX 
NUMBERS  H HOSE  BEAL  AND  IMAGINARY 
PARTS  ARE  STORED  IN  <CPHI>  AND 
<SPHI>  RESPECTIVELY. 

DMTM  -  APPLIES  A  MODIFIED  DIRECT  MELLIN 
TRANSFORM  TO  AN  INPUT  HAVEFORM 
PUTTING  THE  OUTPUT  IN  <XFM>. 

THE  ALGORITHM  IS  BASED  ON  A  FIRST 
BACKHAND  DIFFERENCE. 

SMT  -  APPLIES  A  MODIFIED  MELLIN  TRANSFORM 
BASED  ON  A  SECOND  DIFFERENCE. 

SMT 2  -  APPLIES  A  MODIFIED  MELLIN  TRANSFORM, 
DIFFERENT  THAN  SMT,  BUT  ALSO  BASED 
ON  THE  SECOND  DIFFERENCE. 

CDMT  -  APPLIES  A  MODIFIED  MELLIN  TRANSFORM 

JUST  AS  DMTM  ABOVE,  EXCEPT  THAT  THE 
CENTRAL  DIFFERENCE  IS  USED. 
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SIMP  -  APPLIES  A  MODIFIED  KELLI N  TRANSFORM 
USING  1  BACKHAND  DIFFERENCE  AS  IN 
DMTM,  EXCEPT  THAT  THE  INTEGRATION 
IS  fit  SIMPSON'S  RULE  INSTEAD  OF  THE 
TBAPEZOIDAL  SOLE. 

XAB  -  TAKES  THE  HAGNITUDE  OF  THE  COHPLEX 
TRANSFORM  <XFM>  AND  PUTS  THE 
NANITODE  IN  A  SPECIFIED  VECTOR. 

STON  -  NORMALIZES  A  VECTOR  BY  ITS 

MAGNITUDE  AND  WRITES  IT  TO 
LOGICAL  DEVICE  3  WITH  A  TITLE 
FROM  LOGICAL  DEVICE  <*. 

HOLD  -  LOADS  ONE  VECTOR  INTO  ANOTHER. 

ALTER  -  CHANGES  <WFM>  BT  SCALE  &/OR  SHIFT 
AND  OUTPUTS  TO  A  SPECIFIED  ARRAY 

INTP3  -  A  SECOND  ORDER  SPLINE  INTERPOLATION. 

CFOBH  -  PROVIDES  TWO  CLOSED  FORM  SOLUTIONS 

FOR  VERIFYING  THE  MELLIN  ALGORITHMS. 

TITLE  -  ENTITLES  THE  PLOTS  ON  THE  BASIS  OF 
THE  CALLING  PROGRAM. 


Ill 


BsagliiT  «hin  ~  Ti&v.ttttrrau 


SO  THE  HAIM  PROGRAN  STASTS1 

DIHENSICN  PET (256) ,  BFH  (256)  , STAND  (256)  ,  IF  (5) 
^CCHHON^XFH^256j Ej^CPHI  (256,  128)  , SPHI (256, 128)  ,PI, 

DATA  IP/'  FR' ,*EQ0E»,' NCI  •/ 

PI  3  3.  141592654 


C  H08  HAH I  HAVBP0B9S  ABB  TC  BE  TBANSFOSHED? 

BEAD (2, 10) NUNHFH 
10  POBHAT (14) 


C 

C 

C 

C 

c 

c 


50 


C 

C 


c 

c 

100 

c 


c 

c 


200 


C 

c 

c 

c 


NTHS  IS  THE  HOHBEB  0?  TIHE  SAHPLES  (INCLUDING 
ANY  ZEBO  PILLING)  .  IT  IS  A  PONEB  OF  TNO 
FOB  THE  CONVENIENCE  OP  THE  PPT. 

OPTS  IS  THE  NOMBEB  OF  SAMPLES  INPUT  TO  THE 
HELLIN  TBANSFOBH  BLOCK.  THE  COEFFICIENTS 
ABE  CONFUTES  NON. 

NU  3  8 

8 IMS  3  2** NO 
HPTS  *  NTH S/2 
POBHAT (14) 

NCOEF  *  NTHS 

CALL  CO EF  (NCOEF, HPTS) 

SET  OP  THE  LOOP  FOR  THE  NUHBEB  OF  8AVEF0BHS 
TC  BE  PROCESSED. 

DO  500  IBAVE31 , NUH8FH 

GET  THE  NEXT  INPUT  BAVEFOBH. 

CALL  HAVE  (HFH, NTHS) 

ZEEO  THE  <STAND>  VECTOB  TO  BE  USED  AS  THE 
IHAGINABI  PAST  OF  THE  NTHS  TIHE  SAHPLES. 

DO  100  1*1, NTHS 
STAND  (IT  3  Q.O 
CONTINUE 

CALL  STCB (BFH, NTHS) 

TAKE  THE  FFT. 

CALL  TITLE  (IF) 

CALL  PFT  (HFH, STAND, NTHS ,NG) 

TAKE  THE  HAGNITUDE  AND  PUT  IT  INTO  THE  COHHON 
BAVEFOBH  <HFH>. 

DO  200  1*1, NTHS 

HFH  (I)  3  BFH  (I)**  2  ♦  STAND (I) *• 2 
SQBT  (BFH  (I)  ) 

CONTINUE 

CALL  STOB  (HFH, NTHS) 

TAKE  THE  HEiLLIN  TBANSFOBH  OF  THIS  SPECTBUH 
USING  THE  FIRST  HALF  OF  THE  FFT  SAHPLES, 
OTHEBBISE  KNOWN  AS  HPTS3NTHS/2.  THESE  ABE  THE 
CNLI  UNIQUE  VALUES. 

CALL  DHTH  (8-FM, HPTS,  NCOEF) 

CALL  XABjPBT, NCOEF) 

CALL  STCB  (PBT, NCOEF) 
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NEXT  CEIL  THE  SECOND  ORDER  BULE  DMT 
SUBBOOTINES.  BOTH  COMPOTE  THE  NELLIE 
USING  THE  SECOND  DIEFEBENCE  APPHOXIMATION 
INSTEAD  OF  THE  PIBST  DIEFEBENCE  APPBOXIMATION 
ABOVE.  OTHEBHISE  THE  APPBOACH  IS  THE  SAME. 


CALL  SMI  (HFM, HPTS.NCOEF) 

CALL  XAB (PBT.NCOEF) 

CALL  STCH  (PBT.NCOEF) 

CALL  SHT2  (HFM.MPTS, NCOEF) 

CALL  XAS (PBT.NCOEF) 

CALL  STCH  (PBT, NCOEF) 

CALL  CDBT  HHICH  USES  THE  CENTBAL  EIIFEBENCE 
BOLE  FOB  APFBOXIMATING  THE  TBANSFOBM. 

CALL  CDMT  (HFM.MPTS, NCOEF) 

CALL  XAB (PBT. NCOEF) 

-'#NCOEF) 


CALL  STCH  (PBT, 


CALL  SIMP  HHICH  COMPUTES  THE  MELLIN  TBANSFOBM 
USING  THE  FIBST  DIFFERENCE  ALGORITHM  AND 
SIMPSON *S  BULE  TO  COMPOTE  THE  MODIFIED 
MELLIN  TBANSFOBM. 


CALL  SIMP  (HFM.MPTS, NCOEF) 

CALL  XAB{PBTxNCOBFf 
CALL  STCH  (PBT ,HCQEF) 

END  THE  TBANSFOBM  LOOP.  THE  TBANSFOBM  HAS  BEEN 
OUTPUT  TO  LOGICAL  DEVICE  3  AND  PBEPABED  HIXH 
TITLE  INFORMATION  PROVIDED  BI  LOGICAL  DEVICE  2 
FOB  PLOTTING  HITH  PBOGBAM  MELPLT  FOBTBAN. 

STAX  IN  THE  LOOP  IF  MORE  HA VEFORMS  ABE  AVAILABLE. 
500  CONTINUE 

CALL  CFOBM (PBT, NCOEF) 

CALL  CFOBM (PBT,  NCOEF) 

STOP 

END 
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ft*************************************************: 

*  I HZ  HAVE  SUBROUTINE  PRODUCES  A  WAVEFORM.  IN 

*  BEAD  FROM  THE  LOGICAL  OBIT  2,  NORMALLY  A  SHIP 

*  DATA  PILE.  THE  B  SAMPLES  ABE  BEAD  ABD  THEB  ' 

*  2EBOS  ABE  STUFFED  TO  FILL  THE  256  SAMPLES. 

*  THE  BEQOESTED  WAVEFORM  IS  OUTPUT  IB  THE 

*  COMMON  AfiRAY  WFM (256) . 

***»**«***********tf****************************«**: 


30 

50 


10 

20 


100 


SOBBOUTINS  HAVE  (HFM , NPT S) 

DIMENSION  -  ““  ■" 

COMMON 
1 1XT  ( 1 

DATA 


SIOM  HFM(NPTS) ,ID (5) 

N  XFMj25b#2yCPHI  (256,  128)  ,SPHI  ( 

it/*  BAN'^GE  S*,»  AMPL*,*BS  *,* 


(256,128)  ,PI, 
'/ 


BEAD  (2,  10)  KEY 

- --r 

(IYT 


BEAD  ( 

BEAD 
FOBM 
BEAD (2,  10 


mi! 


BEAD  (2,2o)  1 
CALL  TITLE  ( 
FORMAT  (14) 
FORMAT (F10. 5) 


(I)  ,1*1,10) 

1  #  B) 


BFM(I)  ,1 
ID) 


IF (N  .EQ. 
B  *  N  ♦  1 


NPTS)  BETOBB 


DO  100 
HFM  (I) 
COBTIN 
BETOBB 
END 


I* B, NPTS 

*  o;o 

OB 


114 

-  ^ 


c 

c 

c 

c 

c 

c 

c 

c 

c 


102 


101 

100 


103 


c 

c 

c 

c 

c 

c 


200 


********* ********************************* 

*  SOBBOOTINE  F FT,  GIVES  THE  COMPLEX 

*  SAMPLES  (2**NU=N) ,  SILL  SETUPS  THE 

*  S  COEFFICIENTS  USING  A  EFT. 

****************************************** 


ICO  L»1,NU 

ioi  i*i:s2 


(HI 


SUBBOUTINE  FFT  (XBEAL.XIHAG. N , NO) 
DIMENSION  XBEAL  (N)  ,  XIMAG  (NJ 
H  2*11/2 
B0 1*11 0- 1 
K*0 
DO 

sgiiwMvidii 

C*COS  (A  BG) 

S*SIN  (ABG) 

K1*K+1 
R  1M2*K1 +82 

TBEAL*XBEAL  (K1N2)  «C*XIHAG  (K1S2)  *S 
TIM  AG*XIMAG  (K 1 N2)  *C -XBE  A L  (K 1-N2)  *S 
XBEAL  (K  1N2)  *XBEAl  (K  1)  -TBEAL 
XIMAG  K1N2)  *  XIMAG  iK  li-TIMAG 
XBEAL  K1)  3 XBEAL  (Ki)  +TBEAL 
XIMAG  <K  1)  *XIM AG  (KI)  +TIMAG 
K*K*  1 
K*K*N2 

IF^K.LT.S)  GO  TO  102 

S01*H01— 1 

N2*S2/2 

DO  103  R*1.N 

I»IBITB  (K- 1  cHO)  ♦  1 

IFXI.LE.KIGC  TO  103 

TBEAL*XBEAL(K) 

TIMAG*XIMAG jK) 

XBEAL  K)*XBEAL  (I) 

XIMAG  Kl *XIHAG (I) 

XBEAL  I) *TBEAL 

XIMAG  I1.*IIHAG 

CCBTIBOE 

BZTOBS 

END 


FONCTIOS  IHITB 

FUNCTION  IBITB (J, NU) 

J1*J 

IBITB*0 

DO  200  I* 1 r  NU 

J2-J1/2 

IBITR*IBITB*2+ (J1-2*J2) 
J1*J2 
E  STUBS 
END 
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r 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


************************************************* 

*  THE  COEF  SOBBOUTINE  CONFUTES  THE  NELLIN 

*  COEFFICIENTS  IN  TO  CONHCN  ASSAYS  CPHI  AND 

*  SPHI  THAT  BEPEESENT  THE  REAL  AND  IMAGINARY 

*  PABTS  HESPECTIVELI.  THE  TEBHS  ABE  COHPOTEO 

*  BZ  THE  FOBHOLA: 

*  PHI  (Ir 0)  3  J**S  ,  NHEBE  S  =  A  NOB HAL I ZED 

*  DISCRETE  BAOIAN  FREQUENCY. 

************************************************* 


200 

100 


SOBBOOTINE  COEF  (NC0EF.NPT5) 

CON  HON  XFH  (256, 2)  ,CPHI<256,  12S)  ,  SPHI  (256,  128)  ,PI, 

iiiTnM 


r 


EO  itiu  x  —  i, 

HI  =  FLOAT  (if 
CHEGA  *  2.0 *  PI  * 
00  200  J  »  1.NPTS 
BJ  >  FLOAT  (jf 


BI  /  36. 


CPHI  (I,  J) 

SFHIJI.J) 

CONTINOE 

CONTINOE 

BITOBN 

END 


OS  (ONEGA  * 
a  SIN  (ONEGA  * 


ALOG  (BJ 
ALOG 
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***************************************************** 

I HE  DH1H  SOBBOOTINE  PEBFOBMS  A  DESCBETE  MELLIN 
TBAISFOHH  OH  THE  ABBAY  HFH.  THE  FOBHOLA 
FOB  OHE  HELLIH  FBEQOENCY  FALOE  IS: 

ZFH  (I) *SUH (K=1  TO  NPTS)  (IFM  (1*1) -BFfl(I) ) *K**S 


THE  COHPLEX  COHPONE  NTS  CF  K**S  ABE  COHPOIED  PBIOB 
TC  CALLING  OHTH  AMO  STOBED  IN  THE  COHHOH  ABBAYS 
CPHI  AND  SPHI  AS  BEAL  AND  I MAGI N A BY  PABTS 
BISPECTIYELY.  THE  ALGOBITHH  IS  BASED  ON  THE 
FIBST  DIFFEBENCE.  THE  TBAPEZOID AL  BOLE  IS  OSED 
FOB  THE  INTEGBATION.  THE  COHPLEX  OOTPOT  FOB 
THE  TEA NS FORI  IN  IN  THE  COHHON  ABBAY  <XFH>. 
***************************************************** 


SOBBOOTINE  DHT H (SAHP, NPTS, NCOEF) 
DIMENSION  SAHP (NPTS) , ID  (5) 


DIHBNSION  SAHP (NPTS) .ID  (5) 

CCHHON  XFH (256,2) ,CPHI(256, 128) ,SPHI  (256, 1 
IIXT(IO)  ,IYT  (10)  .KEY 

DATA  ID/' D-HE* , *  LLI N* , *  FBB • , • QOEN' , ' CY  ' 


28),  PI, 


CALL  TITLE  (ID) 

DC  100  1=1, NCOEF 

xfh  (i,i)»o:o 


xfh  (i,ii»o:o 

XFH  fl# 21*0.0 
N 1  -  NPTS  -  1 
CO  200  J*1,N1 
J1  *  J  ♦  1 
XFH  (1,1)  *  CP 
XFHJI.21  *  SP 
CONTINUE 
CCNTINOE 
BETOBN 
END 


CPHI 

SPHI 


(I,J)  * 


(SAHP  (J) 
(SAHP  (J) 


-  SAHP  (J1I)  +  XFH  (1,1) 
SAHP  (Jin  ♦  XFH  (1,2) 


i 


nn  -*  nnnn  no 


******************************************************* 

*  SB?  IS  A  SOBBOUTINE  THAT  COHPOTES  A  NOBEBICAL 

*  APP BO XI BA TICK  TO  THE  BEILIN  TBANSFOBN  AS  DOES 

*  CBSB  ABOVE.  SBT  OSES  A  TBAPAZOIDAL  APPBOXIBAIION 

*  TC  COBPOTE  THE  TBANSFOBN,  BUT  OSES  THE  SANE 

*  COEFFICIENT  BATBIXES  <CPHI>  AND  <SPHI>  CONTAINED 

*  IN  COHHON. 

******************************************************* 


SOBBOOTINE  SBT (SABP,NPTS,NCOEF) 
CIBENSICN  SIBP  (NPTS)  ,ID  (5) 
- - c^I(j5g  - 


CIBENSICN  S-ABP  (NPTS 
COBBON  XFB  (256.21  ,C 
1IXTM0I  ,IYT  (10) -KEY 
CAT!  ID/*  S-nE* , 'LLI 
CALL  TITLE  (ID) 


KEY 

LLIN* #  * 


56,  128)  , SPHI (256, 128) ,PI, 
FBE* ,*Q0EN* , *CY  •/ 


INITIALIZE  THE  INPOT  ABBAY  AND  COBPOTE 
THE  LOOP  CONSTANTS. 

N 1  *  NPTS  -  1 
N2  »  N1  -  1 

SET  OP  THE  TBANSFOBB  LOCP.  THE  OOTEB  LOOP 
SETS  OP  THE  COEFFIECIENTS  WHILE  THE  INNEB 
LOOP  COBPOIES  THE  SOB  WHICH  ABE  THE 
COEFFICIENTS. 


00  200  J  *  1, NCOEF 
XFB  (J,1)  *  Q.O 
XFHjJ-21.  »  Q.O 
DC  100  I  *  1,N2 

10  *  I 

11  *  I  ♦  1 

12  »  I  ♦  2 

DELTA  *  SABWIO)  -  2.*  SAHP  (III  ♦  SABP(I2) 
XFB  (J ,  1 )  *  XFB(J,1)  ♦  DELTA*CPHI  (J,I)*I 

CCMTlfioE  *  XFH’J'2*  *  DELTA*SPHI  (J,I)  *1 

XFB  (0,1)  =  XFB  (J,  11  /  (FLOAT  (J)  *PI/18.) 

XFB  (J ,2)  *  XFB  (J,  2)  /  (FlOiT  (J)  ♦PI/18. ) 


XFB  (J,1) 

XFBJJ-2) 

CONTINUE 

BETOBN 

END 


XFB  (J,  1)  /  SQBT  ( 1+  (FLOAT  (J)  *PI/1 8 . )  **2) 
XFB  (J,2)  /  SQBT  (1+  (FLOAT  (J)  *PI/18.)  **2) 


uuuuuuuuuuu  uu  uuuu 


f 


** ***************************************  **#1(l*#lM»*****,l» 

*  SMT2  IS  A  SUBROUTINE  THAT  COHPUTES  A  NUMERICAL 

*  AP PROXIMATICN  TO  THE  MEX1IN  TBANSFOBM  AS  SOES___ 

*  phtm  ABOVE.  SHT  USES  A  TBAPAZOIDAL  APPROXIMATION 
!  tccoMPUTE  THE  TRANSFORM ,  BUT  USES  THE  SAHE 

*  COEFFICIENT  MATRIXES  <Cpfil>  AND  <SPHI>  CONTAINED 

I*ii***SSJi***********************************#**,M‘<,**# 


SUBBOUTINE_SMT2^SAMP/ NPTSy NCOEF) 


ID  (5 


100 

C 

200 


DIMENSION  SAMP  ... - - —  , 

COMMON  XFM(256. 2j,C&fil  (2 
1IXT  (10)  ,,111  (10;  ,  iCET 
DATA  ID/'  S2-M'  ,  ^  ELLI*  ,  '  N 
CALL  TITLE  (ID) 


,128)  , SPHI  (256,  128)  ,PI, 


FB*  ,  'EQUE'  MCI 


INITIALIZE  THE  INPUT  ABBA?  AND  COMPUTE 
THE  LOOP  CONSTANTS. 

N 1  »  NPTS  -  1 
N2  *  HI  -  1 

SET  UE  THE  TRANSFORM  LOCP.  THE  OUTER  LOOP 
SETS  UP  THE  COEFPIECIENTS  NHILE  THE  INNER 
LOOP  COMPUTES  THE  SOM  MHICH  ARE  THE 
COEFFICIENTS. 


DO  200  J 
XFM  (J,1 
XFH  (J ,2 
DO  100 

10  *  I 

11  *  I  ♦ 

12  »  I  ♦ 


«*  *  ' » 

)  *  0. 

IM: 


1 , NCOEF 
0.0 
0 

N2 


1 

2 


2.*  SAMP  (1 1). 


ifiift~=  =3mi„ 

BiH;H  *  £5  S3:  II  ♦  li:i 

CONTINUE 


♦  SAMP  (12)) 


XFIHJ^I?8*  XPEUJ.1)  /  SORT  (1+  (FLOAT  (J)  *PI/18. 1  ♦* 2) 

mill  -  wmIjUI  '/  s2bt1i*1ploatW*pi/i8.J**2J 

CONTINUE 

RETURN 

END 
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on 


C 

C 

c  ****************************************************** 

C  *  THE  CENT  SUBBOUTINE  FEB F OEMS  A  DESCHETE  MELLIN 
C  *  TBANSFOBM  CM  THE  COMMON  ABB  A I  BFM.  THE  FOBMULA 
C  *  FOB  ONE  M ELLIN  FREQUENCY  VALUE  IS: 

C  *  XFM(I)*SUM(K*1  TO  NETS)  (HFM  (I*  1)-NFM  (1-1)  )  *K**S 

C  *  THE  COMPLEX  COMPONENTS  OF  K**S  ABE  COMPUTED  PBIOB 
C  *  TO  CALLING  CD NT  AND  STOBED  IN  THE  COMMON  ABBAIS 
C  *  CPHI  AND  SPHI  AS  BEAL  AND  IM AGIN AST  PABTS 
C  *  BESPECTIV ELI .  THE  CENTBAL  DIFFEBENCE  IS  USED. 

C  **************************«***************************: 

C 

C 


SUBBOUTINE  CDMT (H, N PIS , NCOEF) 

DIMENSION  H  (NPTS)  ,ID  (S) 

COMMON  XFM  (256 ,2i#CPHI  (256, 128)  ,SPHI  (256,  128)  ,PI, 

hit  ( ioi  ,i yt  (iof  tkki 

DATA  ID/*C-“ME '  ,  '  LLI N*  ,  *  FBE*  ,  '  QUEN*  ,  *  CY  »/ 

CALL  TITLE  (ID) 


C 


100 

C 

200 


SET  UP  DEBXVATIVE  OF  INPUT  VECTOfi  <H (I) > 
IDENTIFIED  HEBE  AS  <G>. 


NG 


NPTS  -  1 


DO  200  J 
XFH  (J,1 
XFH*-  “ 

OMEGA 


1 


0:11  1 8-J 


NCOEF 


FLOAT  (J) 
1,  NG 


*  PI  /  18. 


EC  100  I  » 

IM1  »  I 
IP1  *  I  ♦  2 

DH  >  (H  (I PI)  -  H(IM1))  /  2. 

COMPOTE  THE  J-TH  COEFFICIENTS  BY  THE  SUM. 

13:11 

CONTINUE 


XFH  (J  ,  1 
XFH  * 


8:3 


XFM  (J,1| 
XFH 


DH 

DB 


* 

* 


CPHI 

SPHI 


CCLOB 
XFM  (J,1 

ifh'-  ~ 


8:11 

CONTINUE 


1. 


XFM  (U,1J 
XFM  (J ,  2] 


* 

* 


COLOB 

COLOB 


BETUBN 

END 
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1 


c 

c 

c 

c 

c 

c 

c 


******************************************************* 

*  THIS  SUBBOOTINE  USES  SIMPSON'S  BUIE  TO  APPBOXIMATE 

*  THE  MODIFIED  MELLIN  TBA  NSFOBH.  THE  MODIFICATION 

*  IS  FBEQUENCX  TIMES  THE  FBEQOENCT  DEBIVATIVS  OF 
$  THE  FFx 

******************************************************* 


SUBBOUTIN E  SIMP  (H, NETS, NCOEF) 

EIMENSION  H  (NPTS)  ,  IE  (5) 

COMMON  XFM  (256  ,2}  ,CPHI  (256,  128)  ,SPHI(256,  128)  ,PI, 
1IXT  (10)  ,Ili  (10)  ,  KET  .  . 

DATA  ID/ * I-ME  * , *  LLI N ' , *  FBE • , 'QOEN* , • CT  •/ 

CALL  TITLE  (ID) 

Ml  *  NPTS  -1 
M2  =  NPTS  -  2 


67 

69 


DO  100  J=  1 , NCOEF 
FSIMP  =  2.0 

CMEGA  *  PI  ♦  FLOAT  (J)  /  18. 
COLOB  *  1 . 

XFM  (J,1)  =  Q. 

XFM  (J  ,2)  »  0. 

EC  200  1=1, Ml 
IH1  *  I 
10  =  1  ♦  I 

IF  (FSIMP  .GT. 

FSIMP  *  4. 

GO  TO  69 
FSIMP  »  2. 

CONTINOE 
DELTA  =  H(I0^ 


3.)  GO-  TC  67 


DELTA  =  H  (10)  -  H(IM1) 

XFM  ( J  ,  1 )  =  XFM  (J,  i )  +  FSIMP  *  DELTA  *  CPHI(J,IM1) 
XFM  (J ,2)  =  XFM  (J,  2)  ♦  FSIIJP  *  DELTA  *  SPHI(J,IM1) 


200  COMTIHOE 

100  CONTINUE 
5ZTQBN 
END 
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c 

c  *********«***********•«*********«*********•********) 

C  *  THE  XAB  SUEROUTINE  TAKES  THE  MAGNITUDE  OF  THE  ' 

C  *  COMMON  COMPLEX  ARRAY  XFM(256,2)  AND  PLACES  ' 

C  *  THESE  VALUES  IN  OUTPUT  VECTOR  <XMAG>  FOR  LATER 
C  *  PRINTING. 

C  ***************************************************1 

c 

c 


SUBROUTINE  XAB (XMAG ,NPTS) 

DIMENSION  XHAG(NPTS) 

COMMON  XFM  (256.2)  .CPHI  (256,  128)  , SPHI  (256,  1 28)  , PI, 
1IXTM0)  ,IYT  (10/,KfeY 
DO  100  f=1.NPTS 

*  *QRT  (XFM  (I,  1)  **2+XFM  (1 , 2f **2) 

100  CONTINUE 
EETURN 
END 
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c 

c  ******************************************************* 

C  *  THE  STOii  SUBROUTINE  STORES  THE  INPUT  ABBA  I  PRI(NPTS) 

C  *  INTO  THE  LOGICAL  DEVICE  2  FOB  LATEB  USE  IN  PLOTTING. 

C  *  NPLOl  NUHBEBS  THE  PLOTS  POB  LATEB  IDENTIFICETION . 

C  *  AND  A  PLOT  TITLE  FOB  THE  MELLIN  FREQUENCY  IS  ADDED 

C  *  FOB  CONVENIENCE. 

C  *  THE  HODULUS  OF  THE 

C  *  TBANSFOBH  TAKEN.  SCALED  TO  UNIT  MAGNITUDE,  AND  ' 

C  «  OUTPUT  TO  LOGICAL  UNIT  2. 

C  *  INPUT:  PRT  -  TO  BE  SCALED  TO  1  AND  HBITTEN 

C  *  TO  LOGICAL  UNIT  2. 

C  *  NPLOT  *  THE  NUMBER  OF  THE  PLOT 

C  *  KEY  -  NUMBER  OF  CURVE  THIS  PLOT 

C  ******************************************************* 

C 


CKK 

C 


10 

12 


13 


100 


20 

200 


SUBROUTINE  STON (PfiT , NPTS) 

COHHON^XFMjjlsi), If  f^PHI  (256,  128)  ,  SPHI  (256,  128)  ,PI, 
1IXT  (10)  ,IIT(1C)  ,k6Y 


HBITE  (3,13) NPTS 
iRITE (3, 13V KEY 
““  I  .NE.  1)  GO  TO  12 
AXIS  1ABELS 


IF  (KEY 
HRITE  A 


iRITE  (3, 


0.0 

^a^NPTS 
-  .GT. 


NELOT  » 

BPRT  * 

FORMAT 
DO  100 
IF  (PBT(I) 
CONTINUE 
IF  (BPRT 
DO  200 
PRT  (I)  * 

I  a  FLOAT jl) 
HBITE  (3,20)1 
FORMAT 
CCNTIN 
BETURN 
END 


(IXI  (1 

,1*1,101 

(IYT  (J 

cl  ,1=1.10) 

LOT  ♦ 

1 

BPRT)  BPRT 


.LT.  .00001)  BPBT  a 
1*1, NPTS 

PHT(I)  /  BPRT' 

PET  (I) 

"  ) 


PRT  (I) 

1. 


*9*l{ 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


********************************************** 

*  THE  HOLD  SUBROUTINE  TAKES  THE  INPUT  FILE 

*  <FILIN>  AND  STORES  IT  IN  THE  OUTPUT  FILE 

*  <FILOUT>  FOB  TEflPORARY  STORAGE.  TBE  FILE 

*  <FILIN>  HERRINS  UNCHANGED. 

********************************************** 


SUBROUTINE  HOLD (FILIN, FILOUT, NPTS) 
DIHENSION  FILIN  (NPTS)  ,  FILOUT  (NPTS 


100 


DO  IOC  I  »  1, NPTS 

FILOUT  (I)  *  FILIN  (I) 

CONTINUE 

RETURN 

END 


C 

C 

c 

c 

c 

c 

c 

c 

c 


****************************************************** 

*  THE  ALTER  SUBROUTINE  NHL  ALTER  THE  CORNON  ARRAY 

*  <BFH>  AS  SPECIFIED  BY  TBE  INPUT  VARIABLES 

*  <SCALE>  AND  <SHIFT> .  TBE  ALTERED  HFH  IS  OUTPUT 

*  IN  THE  VECTOR  <ALT>. 

****************************************************** 


SUBROUTINE  ALTER ( ALT, 8F H , SCALE, SHIFT, NPTS) 
DIHENSION  ALT (NPTS) .HFH  (NPTS) .TOLD  (256) .INEH  (256) 
CCHRON  XFR^25o,  2^CEHI  (256,  128)  ,SPHI  (256,128)  ,PI, 


AID  lAJVf 

HIT (10)  ,IYT  (lOf 
DO  100  X* 1 , NPTS 


AW 

FLOAT  (I) 

TOLD  (I)  /  SCALE  ♦ 


SHIFT 


TOLD  (I) 

TNENjli  * 

100  CONTINUE 

CALL  IN TP 3  (HER, TOLD, ALT, THEN, NPTS) 

RETURN 

END 
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nnnnn 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 


******************************************************* 

*  INTP3  IS  A  SECOND  OHDEB  INTERPOLATION  BASED  OB  A 

*  f CLINOMIAL  EQUATING  TO  FIRST  AMD  SECOND  DERIVATIVES 

*  APPBOZI BATED  BI  CENTRAL  DIFFERENCES.  THE  INPUT 

*  VECTOR  <XO>  HAS  OLD  SAHf LES  AT  TINES  IN  <IO>.  THE 

*  MEN  SAMPLE  TIBES  ARE  INPUT  THROUGH  ARRAY  <TN>  AND 

*  THE  CONFUTED  SAMPLES  AT  THESE  TIMES  ARE  OUTPUT  IN 

*  ARE  OUTPUT  IN  THE  VECTOR  <XN>  AND  <HFH>. 
******************************************************* 

SUBROUTINE  INTP3  (XO,TO,  XN.TN.NPTS) 

DIMENSION  X3  (3)  . XO  (NPTS).TO  (NPTSJ.XN  (NPTS)  .  IN  (NPTS) 
COMMON  X?m]256,2) ,CFHI (256, 128)  ,SPHI (256,128) , PI, 
1IXT  (10)  ,ITT  (10)  f&EY 


CHOSE  THE  <TO>  SAMPLE  TIME  CLOSES  TO  THE  HARPED 
TIME  HELD  IN  <TN> 


10 

25 


29 

30 


60 

500 


.AND.  (IN  (I)  .LE.  XPTS)  )  GO  TO  5 


GO  TO  25 


DC  60  1*1 , NPTS 
DT  *  100 

XPTS  =  FLOAT (NPTS) 

HWJH.*84-  '•> 

GO  TO  60 
DO  10  J*1.NPTS 
DTABS  *  ABS  (DT) 

CIST  *  TO(J)  -  TN  (I) 

DIST  *  ABS  (BIST) 

IF  (DTABS  .LI.  DIST) 

DT  *  IN  JI)  -  TO  (J) 

CONTINUE 
JTIHE  »  J  -  1 
JI  *  -1 
DO  30  J*1 ,3 
JI  *  J  ♦  JTIHE-1 

IF  ((JI  .GE.  1)  .AND.  (JI  .LE.  NPTS)) GO  TO  29 
X3(J)  *  0.0 
GO  TO  30 
X3(J)  *  XO(J1) 

CCNTINUB 
TIMM  *  TN  (I). 

TIMQ  *  TO] JTIHE) 

XN  (I)  *  FCN-  (X3,  TIHO  ,TIMN) 

CONTINUE 

CO  500  1*1, NPTS 

CONTINUE 

RETURN 

END 
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r>r>  o  non  nnnn 


tCH  *  A  *  (TN**2)  ♦  3  *  TB  ♦  C 

BETOHB 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


100 

200 

202 

101 


10 

20 


********************************* ***************•******! 

*  THIS  SUBROUTINE  IS  THE  RESULT  0?  A  CLOSED 
«  FORM  CALCULATION  OF  THE  CONTINUOUS  SINC**2 

*  EEING  TAKEN  FROM  THE  TINE  DOMAIN  THROUGH  THE 

*  ENTIRE  FOURIER- MELLIN  PROCESSING.  THE  ALGORITHM 

*  IS  USED  TO  VERIFY  THE  FDM  PROCESS.  THE  OUTPUT 

*  FEATURE  SPACE  SHOULD  BE  IDENTICLE  TO  THAT  PRODUCED 

*  BY  ANY  OF  THE  HELLIN  ALGORITHMS.  FOR  THIS  REASON 

*  TBE  SAMPLE  POINTS  ARE  SYNCHRONIZED  WITH  THOSE  USED 

*  ABOVE.  <NGOEF>  FM  COEFFICIENTS  ARE  USED.  THE 

*  COMMON  VARIABLE  <KEY>  MUST  BE  SET  TO  99  PRIOR  TO 

*  CALLING  TO  GET  THE  SINC**2  OUTPUT  FEATURE  SPACE. 

*  TO  OUTPUT  THE  CLOSED  FORM  SOLUTION  FOR  A  RAMP  IN 

*  THE  FREQUENCY  DOMAIN.  <KEY>  MUST  EQUAL  100. 
******************************************************* 


SUBROUTINE  CFORM (CF ,NCOEF) 

DIMENSION  CF(NCOEF) 

COMMON  ZFM  (256,2) ,CPHI  (256,  128) , SPHI (256, 128) ,PI, 
HIT  (10)  ,IYT  (10)  ,KEY 

IF  (KEY  ,EQ.  100)  GO  TO  200 

IF  (KEY  .ME. 99}  RETURN 

DO  100  M=1,NCOEF 

XM  »  FLOAT  (M)  -  1.0 

HMAG  *  1  ♦  (PI  *  XM  /  18.)**2 

CF  (H)  »  SQRT(HMAG)  /  HMAG 

CONTINUE 

GC  TO  101 

DO  202  I  =  1.NCOEF 
XI  -  FLOAT  (I)  -1.0 
OMEGA  »  PI  •  XI  /  18. 

CFjl£  »  2./SQRT  (4.  ♦ 


COl 


:hue 


OMEGA**2) 


CONTINUE 


BEAD  (2, 20) KEY 

BEAD  (2,  10)  (IXT(I)  ,1  =  1,10) 


FORMAT 


CALL  STOH  (CF,NCOEF) 

RETURN 

END 
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c 

c 

c 

c 

c 

c 


************************************************* 


*  title  titles  the  t  axis  foe  the  plotting 

*  PBOGB1H  ACCORDING  TO  THE  ALGOBITHH  USED  TO 

*  GENERATE  THE  FEATURE  SPACE. 

************************************************* 


SUBROUTINE  TITLE  (ID) 

DIMENSION  ID  (5) 

COM HON  XFH  (256,2) ,CPHI (256,  128)  ,SPHI (256, 128)  ,PI, 
II  XT  (10)  ,  I  XT  (Id/  ,KsX 
DO  fOO  i-1,5 
IYT(I)  »  Ifr(I) 

100  CONTINUE 
RETURN 
END 
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